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ABSTRACT 

We present a description of the adaptive mesh refinement (AMR) implementation of the PLUTO 
code for solving the equations of classical and special relativistic magnetohydrodynamics (MHD and 
RMHD). The current release exploits, in addition to the static grid version of the code, the dis- 
tributed infrastructure of the CHOMBO library for multidimensional parallel computations over 
block-structured, adaptively refined grids. We employ a conservative finite-volume approach where 
primary flow quantities are discretized at the cell-center in a dimensionally unsplit fashion using the 
Corner Transport Upwind (CTU) method. Time stepping relies on a characteristic tracing step where 
piecewise parabolic method (PPM), weighted essentially non-oscillatory (WENO) or slope- limited 
linear interpolation schemes can be handily adopted. A characteristic decomposition-free version of 
the scheme is also illustrated. The solenoidal condition of the magnetic field is enforced by augment- 
ing the equations with a generalized Lagrange multiplier (GLM) providing propagation and damping 
of divergence errors through a mixed hyperbolic/parabolic explicit cleaning step. Among the novel 
features, we describe an extension of the scheme to include non-ideal dissipative processes such as vis- 
cosity, resistivity and anisotropic thermal conduction without operator splitting. Finally, we illustrate 
an efhcient treatment of point-local, potentially stiff source terms over hierarchical nested grids by 
taking advantage of the adaptivity in time. Several multidimensional benchmarks and applications to 
problems of astrophysical relevance assess the potentiality of the AMR version of PLUTO in resolving 
flow features separated by large spatial and temporal disparities. 

Subject headings: hydrodynamics - magnetohydrodynamics (MHD) - methods: numerical - relativity 



1. INTRODUCTION 

Theoretical advances in modern astrophysics have 
largely benefited from computational models and tech- 
niques that have been improved over the past decades. 
In the field of gasdynamics, shock-capturing schemes rep- 
resent the current establishment for reliable numerical 
simulations of high Mach-number, possibly magnetized 
flows in Newtonian or relativistic regimes. As increas- 
ingly more sophisticated methods developed, a number 
of computer codes targeting complex physical aspects to 
various degrees have now become available to the com- 
munity. In the field of magnetohydrodyna mics (MHD), 
examples w orth of no tice are As troBEAR d C unningham 
et al.|| 2009| ), Athena (ISto ne et a l. 2008: Skinner &: Os- 
BATS-R-US (I'i'oth et al. 2011), ECHO (|Dd 
_ 2007), FLASH (i-'ryxcll et al. 20001, NTET 
TTiegler 2008), PLUTO (Mignonc et a 



Zanna et al 
VAN A 



2007 ), 

RAMSEg"(Tcyssie r 20021 [Fromang ct al. 2006TlSl5"A7XC 
([Toth 1996; v an deFHolst et al .||2008|). Some of these 



implementions provide additional capabilities that can 
approach the solution of the equations in the relativis- 
tic regimes: AMRVAC and PLUTO for special rela- 
tivistic hydro, while the ECHO code allows to handle 
general relativistic MHD with a fixed metric. Other 
frameworks were specifically designed for special or gen- 
eral relativistic pur poses, e.g. the RAM code ( Zhang 
& MacFadyen||2006|), HARM (IGammie et al.||2003plSd 
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RAISHIN ( [Mizuno et aLl[2006l ). 

In some circumstances, adequate theoretical modeling 
of astrophysical scenarios may become extremely chal- 
lenging since great disparities in the spatial and tem- 
poral scales may simultaneously arise in the problem of 
interest. In these situations a static grid approach may 
become quite inefficient and, in the most extreme cases, 
the amount of computational time can make the problem 
prohibitive. Typically, such conditions occur when the 
flow dynamics exhibit very localized features that evolve 
on a much shorter scale when compared to the rest of 
the computational domain. To overcome these limita- 
tions, one possibility is to change or adapt the compu- 
tational grid dynamically in space and time so that the 
features of interest can be adequately captured and re- 
solved. Adaptive mesh refinement (AMR) is one such 
technique and can lead, for a certain class of problems, 
to a considerable speed up. Some of the aforementioned 
numerical codes provide AMR implementations through 
a variety of different approaches. Examples worth of 
notice are the pa t ch-bas e d block-structured appro ach of 
Berger & Oligerl (pSll); iBerger & Colellal (11989!) (e.g. 



AS'i'ROBEAR), the tully - octree approach described in 
Dezeeuw &: Powelll ( p93| ; [Khokh lovl (fig gg]) (e.g. RAM- 

([2OOOI) 



Dased octree of MacNcice et al 



SES) or the block- 

(e.g. FLASH) and [Keppens et ahi (i2003j); 
& Keppensl ((2007|) (e.g. BAi'S-R-US, AM 



van der Hoist 

kvAC). 

'i'he present work focuses on the block-structured AMR 
implementation in the PLUTO code and its application 
to computational astrophysical gasdynamics. PLUTO is 
a Godunov-type code providing a fiexible and versa- 
tile modular computational framework for the solution 
of the equations of gasdynamics under different regimes 
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(e.g., classical/relativistic fluid dynamics, Euler/MHD). 
A comprehensive description of the code design and im- 
pl ementation may be fou nd, for the static grid version, 
in Mignone et al. (2007) (paper I henceforth). Recent 



additfons to the code m clude a relativistic ve rsion of the 
HLLD Riemann solver (jMigno ne et al.|[2009[), high-order 
finite difference schemes (Mignone et al.||2U10) and opti 



cally thin ra diative losses with a n on-equilibrium chemi- 
cal network ( Tegileanu et al.|[2008 ). Here we further ex- 
tend the code description and show its performance on 
problems requiring significant usage of adaptively refined 
nested grids. PLUTO takes advantage of the CHOMBO 
librarjrl that provides a distributed infrastructure for 
parallel computations over block-structured adaptively 
refined grids. The choice of block-structured AMR (as 
opposed to octree) is justified by the need of exploiting 
the already implemented modular skeleton introducing 
the minimal amount of modification and, at the same 
time, maximizing code re-usability. 

The current AMR impleme ntation leans o n the Corner- 
Transport-Upwind (CTU, |Colella| [T990t method of 



enceforth) in which 



Mignone &: Tzeferacos] ( |2010| 

a conservative finite volume discretization is adopted 
to evolve zone averages in time. The scheme is di- 
mensionally unsplit, second-order accurate in space and 
time and can be directly applied to relativistic MHD as 
well. Spatial reconstruction can be carried out in primi- 
tive or characteristic variables using high-order interpo- 
lation sc hemes such as the piec ewise parabolic method 



(PPM, IColella fc Woodward| |1984| ), weighted essen- 



tially non-osciliatory (WEiNU) or hnear Total Variation 
Diminishing (TVD) limiting. The divergence-free con- 
straint of magnetic field is e nforced via a mixed h yper- 
bolic/parabolic correction of 'Dedner et al. ([2002) that 
avoids the computational cost associated with an ellip- 
tic cleaning step, and the scrupulous treatment of stag- 
gered fi elds demande d by constrained transport algo- 
rithms (Balsara 2004). As such, this choice provides a 
convenient first step in porting a considerable fraction 
of the static grid implementation to the AMR frame- 
work. Among the novel features, we also show how to 
extend the time-stepping scheme to include dissipative 
terms describing viscous, resistive and thermally con- 
ducting flows. Besides, we propose a novel treatment for 
efficiently computing the time-step in presence of cooling 
and/or reacting flows over hierarchical block-structured 
grids. 

The paper is structured as follows. In Section [2] we 
overview the relevant equations while in Section [3] we de- 
scribe the integration scheme used on the single patch. 
In Section |4] an overview of the block-structured AMR 
strategy as implemented in CHOMBO is given. Sections 
[5] and [6] show the code performance on selected multidi- 
mensional test problems and astrophysical applications 
in classical and relativistic MHD, respectively. Finally, 
in Section [7] we summarize the main results of our work. 

2. RELEVANT EQUATIONS 

The PLUTO code has been designed for the solution 
of nonlinear systems of conservative partial differential 
equations of the mixed hyperbolic/parabolic type. In 
the present context we will focus our attention on the 

https : //seesar . Ibl . gov/anag/chombo/ 



equations of single- fluid magnetohydrodynamics, both in 
the Newtonian (MHD) and special relativistic (RMHD) 



regimes. 



2.1. MHD equations 



We consider a Newtonian fluid with density p, ve- 



locity 



, Vz ) and magnetic induction B 



0, 



{B^, By, Bz) and write the single fluid MHD equations 
as 

dp 

di 

^(pv) 

dt 
d£ 

m + 

dB 

~dt 



V • (pv) 



BB' 



f V-[(f +pt)v-(vB)B] = V-n^-A + pv-g, 
- V X (v X B) = -V X i'qJ) , 



(1) 

where pt = p + B^/2 is the total (thermal+magnetic) 
pressure, £ is the total energy density, g is the gravi- 
tational acceleration term and A accounts for optically 
thin radiative losses or heating. Divergence terms on the 
right-hand side account for dissipative physical processes 
and are described in detail in Section [2.1.11 Proper clo- 
sure is provided by choosing an equation of state (EoS) 
which, for an ideal gas, allows to write the total energy 
density as 



£ = 



P 



1 



1 



r - 1 



-pv 



B^ 



(2) 



with r being the specific heat ratio. Alternatively, by 
adopting a barotropic or an isothermal EoS, the energy 
equation can be discarded and one simply has, respec- 
tivelty, p = p{p) or p — cj.p (where Cs is the constant 
speed of sound) . 

Chemical species and passive scalars are advected with 
the fiuid and are described in terms of their number frac- 
tion Xa where a = 1, • • • , iVjons label the particular ion. 
They obey non-homogeneous transport equations of the 
form 

^^-fV.(pA„v) = p5„, (3) 

where the source term Sa describes the coupling between 
different chemic al elements inside the reaction network 
(see for instance Tegileanu et al.|[2008 1 . 

2.1.1. Non-Ideal Effects 

Non-ideal effects due to dissipative processes are de- 
scribed by the differential operators included on the 
right-hand side of Eq. ([l|. Viscous stresses may be in- 
cluded through the viscosity tensor r defined by 



piy 



Vv + (Vv)^ - -IV • V 



(4) 



where v is the kinematic viscosity and I is the identity 
matrix. Similarly, magnetic resistivity is accounted for 
by prescribing the resistive 77 tensor (diagonal). Dissipa- 
tive terms contribute to the net energy balance through 
the additional flux fl^ appearing on the right hand side 
ofEq. (0: 



T — • J X B . 



(5) 
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where the different terms give the energy flux contribu- 
tions due to, respectively, thermal conductivity, viscous 
stresses and magnetic resistivity. 

The thermal conduction flux Fj, smoothly varies be- 
tween classical and saturated regimes and reads: 



■ class 



(6) 



where q = b(j)pc^ is the magnitude of the saturated flux 
(Cowie & McKee 1977), </) is a parameter of order unity 
accounting for uncertainties in the estimate of g, Ciso is 
the isothermal speed of sound and 



class = K||b( b • Vrj + Kj_ 



VT 



b- vrj 



(7) 



is the classical heat flux with conductivity coefficients 
K|| and K_L along and across the magnetic field lines, re- 
spectively (Orlando et al. 20081. Indeed, the presence 



of a partially ordered magnetic field introduces a large 
anisotropic behavior by channeling the heat flux along 
the field lines while suppressing it in the transverse di- 
rection (here b = B/|B| is a unit vector along the field 
line). We point out that, in the classical limit q ^ oo, 
thermal conduction is described by a purely parabolic 
operator and flux discretization follows standard finite 
difference. In the saturated limit (|VT| — )■ oo), on the 
other hand, the equation becomes hyperbolic and thus 
an upwind discretiza tion of the flux is more appropriate 
( Balsara et al.||2008 ). This is discussed in more detail in 
Appendix |A[ 

2.2. Relativistic MHD equations 

A (special) relativistic extension of the previous equa- 
tions requires the solution of energy-momentum and 
number density conservation. Written in divergence form 
we have 



dt 
dm 

dB 

iH " 
d£ 



-KV-(p7v) = 0, 

- V • [w-i'^^rw - BB - EE] + Vp* = , 

V X (v X B) = , 

V • (m — P7v) = , 



(8) 



where p is the rest-mass density, 7 the Lorentz factor, 
velocities are given in units of the speed of light (c = 
1) and the fluid momentum m accounts for matter and 
electromagnetic terms: m = W7^v + E x B, where E = 
— V X B is the electric field and w is the gas enthalpy. 
Total pressure and energy include thermal and magnetic 
contributions and can be written as 



Pt =P + 



B2-f e2 
2 



£ = wj'^ —p- 



B2 



-P7- (9) 



Finally, the gas enthalpy w is related to p and p via an 
equation of state which can be either the ideal gas law. 



w 



Tp 



or the TM (Taub-Mathews, |Mathews|[l971[ ) equation of 
state 

5 



w 



11 



which provides an anal ytic approximation of the Sy nge 
relativistic perfect gas (Mignone & McKinney 2007). 

A relativistic formulation of the dissipative terms will 
not be presented here and will be discussed elsewhere. 

2.3. General Quasi-Conservative Form 

In the following we shall adopt an orthonormal system 
of coordinates specified by the unit vectors {d is used 
to label the direction, e.g. d — {x, y, z} in Cartesian co- 
ordinates) and conveniently assume that conserved vari- 
ables U — {p, pVjEjB, pXa) - for the MHD equations - 
and U = (p7, m, f , B) - for RMHD - satisfy the following 
hyperbolic/parabolic partial differential equations 



dlJ 

~dt 



V • F = V • n + S 



p ' 



(12) 



where F and □ are, respectively, the hyperbolic and 
parabolic flux tensors. The source term Sp is a point- 
local source term which accounts for body forces (such 
as gravity), cooling, chemical reactions an d the source 
term for the scalar multiplier (see Eq. 14 below). We 



note that equations containing curl or gradient opera- 
tors can always be cast in this form by suitable vector 
identities. For instance, the projection of V x E in the 
coordinate direction given by the unit vector can be 
re- written as 

(V X E) • = V • (E X Gd) + E • (V X e<j) , (13) 

where the second term on the right hand side sh ould 
be included as an additional source term in Eq. (12 1 



(10) 



whenever different from zero (e.g. in cylindrical geome- 
try). Similarly one can re- write the gradient operator as 
Vp = V • {\p). 

Several algorithms employed in PLUTO are best 
implemented in terms of primitive variables, V = 
(p, v,B,p). In the following we shall assume a one-to- 
one mapping between the two sets of variables, provided 
by appropriate conversion functions, that is V = V(U) 
and U = U(V). 

3. SINGLE PATCH NUMERICAL INTEGRATION 

PLUTO approaches the solution of the previous sets 
of equations using either finite-volume (FV) or finite- 
difference (FD) methods both sharing a flux- conservative 
discretization where volume averages (for the former) or 
point values (for the latter) of the conserved quantities 
are advanced in time. The implementation is based on 
the well-established framework of Godunov type, shock- 
capturing schemes where an upwind strategy (usually 
a Riemann solver) is employed to compute fluxes at 
zone faces. For the present purposes, we shall focus on 
the FV approach where volume-averaged primary flow 
quantities (e.g. density, momentum and energy) retain 
a zone-centered discretization. However, depending on 
the strategy chosen to control the solenoidal constraint, 
the magnetic field can evolve either as a cell-average or 
as a face-average quantity (using Stoke's theorem). As 
described in paper I, both approaches are possible in 
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PLUTO by choosing between Powell's eight wave for- 
mulation or the constrained transport (CT) method, re- 
spectively. 

A third, cell-centered approach based on the general- 
ized L agrange multiplier (GLM) formulation of Dedner 
et al. |2002) has recently been introduced in PLU'i'U and 
a thorough discussion as well as a direct comparison with 
CT schemes can be found in the recent work by MT. The 
GLM formulation easily builds in the context of MHD 
and RMHD equations by introducing an additional scalar 
field ijj which couples the divergence constraint to Fara- 
day's law according to 



dB 

~dt 

I ~dt 



- V X (v X B) + V?/; 
+ c2 V • B 



= 0, 



(14) 



where Ch is the (constant) speed at which divergence er- 
rors are propagated while Cp is a constant controlling the 
rate at which monopoles are damped. The remaining 
equations are not changed and the conservative charac- 
ter is not lost. Owing to its ease of implementation, we 
adopt the GLM formulation as a convenient choice in the 
development of a robust AMR framework for Newtonian 
and relativistic MHD flows. 



3.1. Fully Unsplit Time Stepping 

The system of conservation laws is adv anced in. time 
using the Corner- Transport-Upwind (CTU, Colella 1990 1 
method recently described by MT. Here we outline the 
algorithm in a more concise manner and extend its ap- 
plicability in presence of parabolic (diffusion) terms and 
higher order reconstruction. Although the algorithms il- 
lustrated here are explicit in time, we will also consider 
the description of more sophisticated and effective ap- 
proaches for the treatment of parabolic operators in a 
forthcoming paper. 

We shall assume hereafter an equally-spaced grid with 
computational cells centered in {xi,yj,Zk) having size 
Ax X Ay X Az. For the sake of exposition, we omit 
the integer subscripts i, j or k when referring to the cell 
center and only keep the half-increment notation in de 



noting face values, e.g., F^^-i- = Fy^^ j_^_i j. when d = y. 
Following 'Mignone et all (|2007|) , we use AV and Ad,± 



to denote, respectively, the cell volume and areas of the 
lower and upper interfaces orthogonal to e^. 

An expl icit second-order accurate discretization of 
Eqns. (12), based on a time-centered flux computation, 
reads 



U" + Ai" 



E 



"PA 



(15) 



where U is the volume-averaged array of conserved val- 
ues inside the cell {i,j,k), At" is the explicit time step 
whereas CH,d and Cp,d are the increment operators cor- 
responding to the hyperbolic and parabolic flux terms. 



respectively: 



r"^2 _ _ 
'-H,d ^ 



n+i 
d.-' 



AV, 



'-'pa'' ^ ' 



n+i 



AVd 



(16) 



(17) 



In the previous expression 'Pd.± and Il^.i are, respec- 
tively, right (-I-) and left (— ) face- and time-centered ap- 
proximations to the hyperbolic and parabolic flux com- 
ponents in the direction of ed- The source term rep- 
resents the directional contribution to the total source 
vector ^ Sd = Sbody + Sgoo including body forces and 
geometrical terms implicitly arising when differentiating 
the tensor flux on a curvilinear grid. Cooling, chem- 
ical reaction terms and the source term in Eq. |14| 
(namely (c^/Cp)V') are treated separately in an operator- 
split fashion. 

The computation of 'Pd,± requires solving, at cell in- 
terfaces, a Riemann problem between time-centered ad- 
jacent discontinuous states, i.e. 



p"+5 


= n 






= n 






^ n 





'n+i 



n+i 
k,+ 



' ^i+l.- 



(18) 



where 7^(-, •) is the numerical flux resulting from the solu- 
tion of the Riemann problem. With PLUTO , different 
Riemann solvers may be chosen at runtime depending 
on the selected physical module (see Paper I): Rusanov 
(Lax-Friedrichs) , HLL, HLLC and HLLD are common to 
both classical and relativistic MHD modules while the 
Roe solver is available for hydro and MHD. The input 
states in the Eq. ( 18 1 are obtained by first carrying an 



evolution step in the normal direction, followed by cor- 
recting the resulting values with a transverse flux gra- 
dient. This yields the corner-coupled states which, in 
absence of parabolic terms (i.e. when fl = 0), are con- 
structed exactly as illustrated by MT. For this reason, 
we will not repeat it here. 

When n 0, on the other hand, we adopt a slightly 
different formulation that does not require any change 
in the computational stencil. At constant x-faces, for 
instance, we modify the corner-coupled states to 



Tt"+5 _ 



u* 



At""' 



d^x d 



(19) 



where, using Godunov's flrst-order method, we compute 
for example 



■ y,j+i 



7^(u^,u^+l), 



(20) 



i.e., by solving a Riemann problem between cell-centered 
states. States at constant y and z-faces are constructed 
in a similar manner. The normal predictors U^- can be 
obtained either in ch aract eristic or primitive variables 
as outlined in Section 3.2 and Section |3.3[ respectively. 
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Parabohc (dissipative) terms are discretized in a flux- 
conservative form using standard central finite differ- 
ence approximations to the derivatives. For any term 
in the form 11 — g{'U)dxf{V), for instance, we evaluate 



the right interface flux appearing in Eq. (17) with the 
second-order accurate expression 



n 



U„; + U, 



i+1 



/(u. 



l+l) 



(21) 



and simila rly for the other d irections (a similar approach 
is used by Toth et aL]|2008 1 . We take the solution avail- 
able at the cell center at t = in Eq. ( 19 1 and at the 
half time step n -|- ^ in Eq. (15). For this latter up- 
date, time- and cell-centered conserved quantities may 
be readily obtained as 



U 



n+i 



(22) 



The algorithm requires a total of 6 solutions to the Rie- 
mann problems per zone per step and it is stable under 
the Courant-Friedrichs-Levy (CFL) condition Ca < 1 (in 
2D) or Ca < 1/2 (in 3D), where 



Ca = At 



n+l 



max 

ijk 



max 

d 



\ max 


Axd 



E 



2pmax 



(23) 



is the CFL number. A™'''' and P^^x ^^^^ ^j^g (local) largest 
signal speed and diffusion coefficient in the di rect ion 
given by c? = {x,y,z}, respectively. Equation (23) is 



used to retrieve the time step At"+^ for the next time 
level if no cooling or reaction terms are present. Other- 
wise we further limit the time step so that the relative 
change in pressure and chemical species remains below a 
certain threshold e^: 



At" 



At" 



ecAt" 
max{\5p/p\,\dX^\) 



(24) 



where 6p/p and SX,^ are the maximum fractional v aria- 
tions during the source step ( Te§ileanu et al.pOOS ). We 
note that the time step limitation given by equation ( 24 1 



does not depend on the mesh size and can be estimated 
on unrefined cells only. This allows to take full ad vant age 
of the adaptivity in time as explained in Section 4.2 



3.2. Normal predictors in characteristic variables 

The computation of the normal predictor states can be 
carried out in characteristic variables by projecting the 
vector V of primitive variables onto the left eigenvectors 
If = l''(Vi) of the primitive system. Specializing to the 
x direction: 



If • V 



i+l 



I 



-s, 



(25) 



where k = 1, • • • , A'^vave labels the characteristic wave 
with speed A** and the projection extends to all neighbor- 
ing zones required by the interpolation stencil of width 
25" + 1. The employment of characteristic fields rather 
than primitive variables requires the additional compu- 
tational cost associated with the full spectral decompo- 
sition of the primitive form of the equations. Neverthe- 
less, it has shown to produce better-behaved solutions 



for highly nonlinear problems, notably for higher order 
methods. 

For each zone i and characteristic field k, we first inter- 
polate wfi to obtain wfj., that is, the rightmost (-I-) and 
leftmost (— ) interface values from within the cell. The in- 
terpolation can be carried out using either fourth-, third- 
or second-order piecewise reconstruction as outlined later 
in this section. 

Extrapolation in time to t" + At"/2 is then carried out 
by using an upwind selection rule that discards waves 
not reaching a given interface in At/2. The result of 
this construction, omitting the k index for the sake of 
exposition, reads 



.,rcf 



[+13+ {w,^+ ~ ^ [Sw, + 6^w, (3 - 2iy)] - , 

(26) 

,f_+l3- {w,,- - \ {&w, ~ S^w, (3 + 21^)] - <!} , 

(27) 

where v = XAt/Ax is the Courant number of the K-th 
wave, /3± = (1 ± sign(i^))/2, whereas Swi and S'^Wi are 
defined by 



Swi = 



S'^Wi = Wi_+ — 2wi + Wi,^ . (28) 



is somewhat arbi- 



The choice of the reference state 
trary and one can simply set wl^^ = Wifi (Rider et al 



20071 which has been found to work well for flows con 



tammg strong discontinui ties. Alternatively, one can use 
the original prescription ( [Colella fc Woodward||1984 ) 



,.rcf 



w. 



rcf 



[Swi + S'^Wi (3 - 
[Swi - 5^Wi (3 - 



2Vrr 



0] 



(29) 
(30) 



max(0, maxK;(z^K;)) and 



where 

min(0, min„(z^K;)) are chosen so as to minimize the size 



of the term susceptible to characteristic limiting ( Colella 
& W oodward! 1 1984) [Miller fc CoTeIIa| [20021 IMignone et 



|al.||2005p . However we note that, in presence of smooth 
hows, both choices may reduce the formal second order 
accuracy of the scheme since, in the limit of small At, 
contributions carried by waves not reaching a zone edge 
are not include d when reconstruc ting the corresponding 
interface value ( Stone et al.[[2008 ). In these situations, a 
better choice is to construct the normal predictors with- 
out introducing a specific referen ce st ate or . eq uivalentlv. 



by assigning u'J'^^ = Wi^± in Eq. \2Q\ and ([27 1 



The time-centered interface values obtained in charac 



teristic variables through Eq. (26) and (27) are finally 



used as coefficients in the right-eigenvector expansion to 
recover the primitive variables: 



At" 



ABx 
' Ax 



'f'^' Ax 



(31) 

where is the right eigenvector associated to the K-th 
wave, Sg is a source term accounting for body forces and 
geometrical factors while Sb^ and arise from calcu- 
lating the interface states in primitive variables rather 
than conservative ones and are essential for the accuracy 
of the scheme in multiple dimensions. See MT for a de- 
tailed discussion on the implementation of these terms. 



6 



Mignone et al. 



As mentioned, the construction of the left and right 
interface values can be carried out using different inter- 
polation techniques. Although some of the available op- 
tions have already be en presented in the original paper 
([Mignone et al. 2007), here we briefly outline the imple- 
mentation details tor three selected schemes providing 
(respectively) fourth-, third- and second-order spatially 
accurate interface values in the limit of vanishing time 
step. Throughout this section we will make frequent us- 
age of the undivided differences of characteristic variables 
p5|) such as 



Aw, 



1^4, + 1 



Wj,0 ■ 



Al 



Win 



(32) 



for the i-th zone. 



Piecewise Parabolic Method — The origi nal PPM recon- 
struction by Colella & Woodward ( 1984 ) (see also Miller' 
[fc ColcUa 2002| r can be directly apphed to characteristic 
variables givmg the following fourth-order limited inter- 
face values: 

w^,± = ^ T Q , (33) 

where slope limiting is used to ensure that 'Wi^± are 
bounded between Wi and Wi±i: 



A 



Wi,o = mm 



Aw^ +1 -I- Aw, 
2 



- , 2mm ( Aw, _ i , Aw. 



,h) 



and 



mm(a, b) 



sign(a) + sign(b) 



(34) 



.(|a|,|6|), (35) 



is the MinMod function. The original interface values 
defined by Eq. ( 33 1 must then be corrected to avoid the 



appearance of local extrema. By defining S± — Wi^± — 
Wifi, we further apply the following parabolic limiter 



S+ = 



if (5+<5_ > 

-26^ if \6±\>2\S^\, 



(36) 



where the first condition flattens the distribution when 
Wi^Q is a local maximum or minimum whereas the second 
condition prevents the appearance of an extremum in 
the parabolic profile. Once limiting has been applied, 
the final interface values are obtained as 



(37) 



In some circumstances, we hav e fou nd that further appli- 
cation of the parabolic limiter ( 36 1 to primitive variables 
may reduce oscillations. 

Third-order improved WENO — As an alternative to 
the popular TVD limiters, the third-order improved 
weighted essentia lly non-oscillatory (WENO) re construc- 
tion proposed by Yamaleev fc Carpenter ( 2009 1 (see also 
Mignone et al. 2010) may be used. I'he interpolation still 
employs a three-point stencil but provides a piecewise 
parabolic profile that preserves the accuracy at smooth 
extrema, thus avoiding the well known clipping of classi- 
cal second-order TVD limiters. Left and right states are 
recovered by a convex combination of 2" order inter- 
polants into a weighted average of order 3. The nonlin- 
ear weights are adjusted by the local smoothness of the 



solution so that essentially zero weights are given to non 
smooth stencils while optimal weights are prescribed in 
smooth regions. In compact notation: 



W' 



Wi, 



where 



W^n 



Wi.o 



At 



2a+ -|- a_ 
a-Aw,_i + ia+Am, +1 



(38) 



a± = 1 



2a. 



a+ 



Al 



Al 



(39) 



As one can see, an attractive feature of WENO recon- 
struction consists in completely avoiding the usage of 
conditional statements. The improved WENO scheme 
has enhanced accu racy with respect t o the traditional 



3rd order scheme of Jiang & Shu 



1996 ) in regions where 



the solution is smooth and provides oscillation-free pro 
files near strong discontinuities . 



Linear reconstruction - 
ing is provided by 



Second-order traditional limit- 



Wi+ = W,o ± 



Awi^o 



(40) 



where Awi n is a standard limiter function such as the 



monotonized-central (MC) limiter (Eq. 34). O ther, less 



steep forms of limiting are the harmonic mean ( van Leer 
T974| 



Awi 



2Awi_^_iAw, 
Aw, ,1+ An 



if Al 



> 0, 



otherwise ; 



the Van Albada limiter ( van Albada et aL]|1982 ): 



i(Awi + i +Aw,^^i) 



Aw 



i,a 



Aw^ , 1 + Aw'^ 



(41) 



if Aw, 







otherwise 
(42) 



or the MinMod limiter, Eq. ([351. 

Linear reconstruction may also be locally employed in 
place of a higher order method whenever a strong shock 
is detected. This is part of a built-in hybrid mechanism 
that selectively identifies zones within a strong shock in 
order to introduce additional dissipation by simultane- 
ously switching to the HLL Riemann solver. Even if the 
occurrence of such situations is usually limited to very 
few grid zones, this fail-safe mechanism does not sac- 
rifice the second-order accuracy of the scheme and has 
been found to noticeably improve the robustness of the 
algorithm avoiding the occurrence of unphysical states. 
This is described in more detail in AppendixiBj However, 
in order to assess the robustness and limits of applicabil- 
ity of the present algorithm, it will not be employed for 
the test problems presented here unless otherwise stated. 

3.3. Normal predictors in primitive variables 
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The normal predictor states can also be directly com- 
puted in primitive variables using a simpler formulation 
that avoids the characteristic projection step. In this 
case, one-dimensional left and right states are obtained 

at t"+5 



usmg 



V*, = V 



v(u*)-vr± 



AV,; 



(43) 



where /3+ = (1-1- sgn(Amax))/2 and /3_ = (1 - 
sgn(Ai„in))/2 may be used to introduce a weak form 
of upwind limiting, in a similar fashion to Section |3.2[ 
Time-centered conservative variables U* follow from a 
simple conservative MUSCL-Hancock step: 



U* = U"- 



(44) 

where Ax ± and AV are the area and volume elements 
and V" _i_ are obtained using linear reconstruction (Eq. 
40 ) of the primitive variables. This approach offers ease 
of implementation over the characteristic tracing step as 
it does not require the eigenvector decomposition of the 
equations nor their primitive form. Notice that, since 
this step is performed in conservative variables, the mul- 
tidimensional terms proportional to dB^/dx do not need 
to be included. 

4. AMR STRATEGY - CHOMBO 
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Figure 1. Two-dimensional example of a three-level AMR hi- 
erarchy, with the base level (£ = 0) covering the entire computa- 
tional domain. Solid lines are representative of the level resolution. 
Dashed lines contour the ghost zones of two patches of level 1 = 1. 
Colors indicate different filling methods: physical outer boundaries 
(red), boundaries filled by exchanging values with adjacent patches 
on the same level (blue), boundaries filled by interpolating from the 
next coarser level (yellow). 



The support for Adaptive Mesh Refinement (AMR) 
calculations in PLUTO is provided by the CHOMBO 
library. CHOMBO is a software package aimed at pro- 
viding a distributed infrastructure for serial and paral- 
lel calculations over block-structured, adaptively refined 
grids in multiple dimensions. It is written in a combi- 
nation of C-I--I- and Fortran77 with MPI and is devel- 
oped and distributed by the Applied Numerical Algo- 
rithms Group of Lawrence Berkeley National Laboratory 
(https://seesar.lbl.gov/anag/chombo/). 



In the block-structured AMR approach, once the so- 
lution has been computed over a rectangular grid which 
discretizes the entire computational domain, it is possible 
to identify the cells which require additional resolution 
and cover them with a set of rectangular grids (also re- 
ferred to as blocks or patches), chara cterized by a finer 
mesh spac ing. CHOMBO follows the|Berger fc Rigout-| 



SOS 



( 1995 1 strategy to determine the most efficient patch 
layout to cover the cells that have been tagged for refine- 
ment. This process can be repeated recursively to define 



the solution over a hierarchy of ^ = 0, . 



levels of 



refinement whose spatial resolutions satisfy the relation 
Acc^ = Ax^"''^, where the integer is the refinement 
ratio between level ^ and level £+\. A level of refinement 
is composed by a union of rectangular grids which has 
to be: disjointed, i.e. two blocks of the same level can 
be adjacent without overlapping; properly nested, i.e. a 
cell of level £ cannot be only partially covered by cells 
of level i +1 and cells of level £ + 1 must be separated 
from cells of level ^ — 1 at least by a row of cells of level 
A simple example of a bidimensional adaptive grid 
distributed over a hierarchy of three levels of refinement 
is depicted in Fig. [T] 



Following the notation of Pember et al. ( 1996 1, a global 
mapping is employed on all levels: m three dimensions, 
cells on level L are identified with global indexes i,j,k 
(0 < i < Nl < j < N^, < k < Ni, with 7Vf_2,3 being 
the equivalent global resolution of level £ in the three 
directions). Correspondingly, the cell i,j,k of level £ is 
covered by (r^)^ cells of level £ + 1 identified by global 
indexes I, m, n satisfying the conditions r^i < I < r^{i -\- 
1) - 1, r^j <m< r\j + l)-l,r^k<n< r^(k + 
Taking direction 1 as an example, the expressions x\ _ = 
iAx^, x{ = {i + 1/2) Aa:^, + = {i + l)^x{, define the 
physical coordinates of the left edge, center and right 
edge of a cell respectively. 



Level 2 



Level 1 



Level 





to+dt/4 



to+dt/2 



to+3dt/4 



t„+dt 



Figure 2. Schematic representation of the time evolution of an 
AMR hierarchy composed by three levels with a refinement ratio 
r = 2. The length of the horizontal black arrows is proportional 
to the timestep size At*. Curved vertical arrows indicate interlevel 
communications. Red arrows represent fine-to-coarse communica- 
tions between s ync hronized adjacent le vels, in cluding conservative 
averaging (Eq. |47| and refluxing (Eq. |48[ |49[l . Green arrows rep- 
resent coarse-to-fa ne c ommunications, mcludmg the conservative 
interpolation (Eq. |45| needed to fill ghost zones and to define the 
solution on newly generated cells of the finer level. 

If the adaptive grid is employed to evolve time- 
dependent hyperbolic partial differential equations, the 
CFL stability condition allows to apply refi nement in 
time a. s well as in space, as first propos ed by Berger fc 
Oliger (1984) and further developed in jBerger &: Colella 
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Table 1 

Systems of coordinates adopted 
in PLUTO-CHOMBO 





Cartesian 


Cylindrical 




X 


r 




y 


z 




z 


1 




X 


r2/2 


y2 


y 


z 


y-i 


z 


1 




AyAz 


r^Az 




AzAx 


rAr 




AxAy 


1 


V 


AxAyAz 


rArAz 



( 1989 ). In fact, each level I advances in time with a time- 
step Ai^ = /S.t'^^^ / r^^^ which is r^^^ times smaller than 
the time-step of the next coarser level £ — 1. Starting 
the integration at the same instant, two adjacent levels 
synchronize every timesteps, as schematically illus- 
trated in Fig. [2^ for a refinement ratio = 2. Even 
though in the following discussion we will assume that 
level -t- 1 completes timesteps to synchronize with 
level CHOMBO allows the finer level ^ -f 1 to advance 
with smaller substeps if the stability condition requires 
it. Anyway, the additional substeps must guarantee that 
levels £ and are synchronized again at time t^ + At^. 

The time evolution of single patches is handled by 
PLUTO , as illustrated in Section [3j Before starting the 
time evolution, the ghost cells surrounding the patches 
must be filled according to one of these three possibili- 
ties: (1) assigning "physical" boundary conditions to the 
ghost cells which lie outside the computational domain 
(e.g. the red area in Fig. [T]); (2) exchanging boundary 
conditions with the adjacent patches of the same level 
(e.g. the blue area in Fig. fl|; (3) interpolating values 
from the parent coarser levelior ghost cells which cover 
cells of a coarser patch (e.g. the yellow area in Fig. [I]). 



Level i+A 




Figure 3. Illustrative bidimensional example of prolongation and 
restriction operations between two levels with a refinement ratio 
= 2. Cells on level i + 1 can be filled by linearly interpolating 
the coarse values on level £ (empty cir cles ) at cross-marked points 
(prolongation, green dashed lines, Eq. |45| . Cells on level £ can be 
filled by averaging down values from leve l ^ + 1 in a conservative 
way (restriction, red dotted lines, Eq. |47[l. 



main program 

initialize(0 l^,^) 

for g = 0, iVsieps 

advance(0, qAt". At") 
endfor 

end main program 



Assign initial conditions on entire level hierarchy. 



Advance base level for lYsteps timesteps. 



procedure advance(f , A^^) 

if (need_to_regrid) 
tag_and_regrid(^) 
endif 



Tag cells for refinement (Eq. 50) and generate new 
grids on levels {! + 1 ^'max- Fill new cells inter- 
polating from coarse levels (Eq. 45). 



fill.boundaries(U', U'l*) 



evolve(U'. A;') 

if (e < fi„a.x) then 
Af'+i = At'/r'' 
for g = 0, . . . , r'-l 
ad«ance(f +l,t' + qAt'+'^ 
endfor 

average_down(U^+\ U^) 
reflux(U',iF^) 
endif 
end procedure 

Figure 4. Pseudocode 
block structured AMR. 



Fill ghost zones on level (: physical boundaries, 
exchange values with adjacent patches, interpolate 
from level f - 1 (Eq. 45 and 45). 

Advance in time level ( from t'- to t^- + At^ (Eq. 
15) and store fluxes at coarse-fine boundaries for 
conservative refluxing. 



At^'^^) Advance in time level f+1 for timesteps. 



Average level [ + I down to level ( (Eq. 47). 

Add conservative flux correction (Eq. 48 and 49) to 
level I 



for the recursive level integration in the 



As schematically illustrated in Fig. [3] in two dimen- 
sions, the coarse-to-fine prolongation needed in case (3) 
(green dashed lines) is based on a piecewise linear in- 
terpolation at points marked by crosses using the linear 
slopes computed from the stirrounding coarse cells. In 
three dimensions the interpolant has the general form: 



U 



U: 



3 

E 

d=l 



i+1 



V, 



V5 



(45) 



where is the volume coordinate of cell centers of level 
£ in direction d and _|_ is its value on the right and 
left faces of the cell respectively (see Table [T] for def- 
initions) . The linear slopes A^Uf ^ are calculated as 
central differences, except for cells touching the domain 
boundary, where one-sided differ enc es are employed. The 
monotonized-central limiter (Eq 34 ) is applied to the lin- 



ear slopes so that no new local extrema are introduced. 

Notice that, since two contiguous levels are not always 
synchronized (see Fig. [2]), coarse values at an intermedi- 
ate time are needed to prolong the solution from a coarse 
level to the ghost zones of a finer level. Coarse values of 
level £ are therefore linearly interpolated between time 

an d ti me + At^ and the piecewise linear interpolant 
Eq. (|45|) is applied to the coarse solution 



Uf* = (1 - a)VlJt') + aVl^,{t' + At') , (46) 
where a = (i^+^ — t')/ At'. This requires that, everytime 
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level £ and level i+1 are synchronized, a timestep on level 
£ must be completed before starting the time integration 
of level £+1. Therefore, the time evolution of the entire 
level hierarchy is performed recursively from level £ = 
up to ^ = ^niaxj ELS Schematically illustrated by the 
pseudo-code in Fig. [4] 

When two adjacent levels are synchronized, some cor- 
rections to the solutions are performed to enforce the 
conservation condition on the entire level hierarchy. To 
maintain consistency between levels, the solution on the 
finer level £ + 1 is restricted to the lower level £ by aver- 
aging down the finer solution in a conservative way (red 
dotted lines in Fig. [3| 



after the time integration of level £ and 
completed: 



V 



1 has been 



(49) 



Finally, when levels from £ up to ^max are synchronized, 
it is possible to tag the cells which need refinement and 
generate a new hierarchy of grids on levels from ^ -I- 1 up 
to £max which covers the tags at each level. Whenever 
new cells are created on level -I- 1 it is possible to fill 
them interpolating from level £ according to Eq. ( 45 1 . It 



U: 



E 



Ir-^O + l)- 

E 

rn—r^j 



lr'^{k+l)- 

E 

n—r^k 



u 



i+1 



V: 



where Vlj ,, 



(47) 

is the volume of cell i,j,k of level £. 
Moreover, the flux through an edge which is shared 
between a cell of level £ and (r^)^ cells of level £ + 1 
must be corrected to maintain the conservative form of 
the equations. For example if the cell i,j,k on level £ 
shares its left boundary with (r^)^ cells of level £+1 (see 
Fig. [5]) , the flux calculated during the coarse integration 
must oe replaced with the average in time and space of 
the fluxes crossing the (r^)^ faces of the finer level cells. 
In this particular example, the flux correction is defined 



is important to notice that this interpolant preserves the 
conservative properties of the solution. 

4.1. Refinement Criteria 

In PLUTO-CHOMBO zones are tagged for refinement 
whenever a prescribed function x(U) of the conserved 
variables and of its derivatives exceeds a prescribed 
threshold, i.e., x(U) > Xr- Generally speaking, the re- 
finement criterion may be problem-dependent thus re- 
quiring the user to provide an appropriate definition of 
x(U). The default choice ado pts a criterion based on the 
second derivative error norm (L6hner||1987), where 



X(U) 



(50) 



where a = a{\J) is a function of the conserved variables, 
± 1 are the undivided forward and backward differ- 



^ E K^+^d^rt%,+ ' '^^'^^^ direction d, e.g., A^^icr = ±(cr,±i 



(48) 

where the index q sums over the timsteps of level £ -\- \, 
while m and n are the indexes transverse to direction d. 
The flux correction is added to the solution on level £ 



o-j). 

The last term appearing in the denominator, ad. rRU pre- 
vents regio ns of small ripples from being refined ( Fryxell| 
et al.||2000[ ) and it is defined by 



|cr,+l| + 2|cri| + |cr,_i| 



(51) 




Similar expressions hold when d — y or d — z. In the 
computations reported in this paper we use e — 0.01 as 
the default value. 

4.2. Time Step Limitation of Point-Local Source 
Terms 

In the usual AMR strategy, grids belonging to level £ 
are advanced in time by a sequence of steps with typical 
size 



mm 



i2y 



(52) 



where we assume, for simplicity, a grid jump of 2. Here 
Ai^jjjj is chosen by collecting and re-scaling to the base 
grid the time steps from all levels available from the pre- 
vious integration step: 



At" 



min r(2)^At^'"] , 



(53) 



Figure 5. Schematic visualization of the refluxing operation 
needed at fine-coarse interfaces to preserve the conservative prop- 
erties of the solution. One cell on the coarser level £ and the cells 
of level £ + 1 adjacent on one side are represented, assuming a 
refinement ratio = 2. Whenever level £ and £ + 1 are synchro- 
nized, the coarse flux (red arrow) must be replaced by the spatial 
and temporal average of t he finer fluxes crossing the fine-coarse 
interface (blue ar row s, Eq. |48| l and the solution must be corrected 
accordingly (Eq. |49[|. 



where At '" is computed using (24). However, this pro 



cedure may become inefficient in presence of source terms 
whose time scale does not depend on the grid size. As 
an illustrative example, consider a strong radiative shock 
propagating through a static cold medium. In the opti- 
cally thin limit, radiative losses are assumed to be lo- 
cal functions of the state vector, but they do not in- 
volve spatial derivatives. If the fastest time scale in the 
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problem is dictated by the cooling process, the time step 
should then become approximately the same on all lev- 
els, At « At^j^j « At°^j, regardless of the mesh size. 
However from the previous equations, one can see that 
finer levels with £ > will advance with a time step (2)^ 
sma ller than required by the single grid estimate. Eq. 
(52) is nevertheless essential for proper synchronization 



between nested levels. 

Simple considerations show that this deficiency may 
be cured by treating split and leaf cells differently. Split 
zones in a given level £ are, in fact, overwritten dur- 
ing the projection step using the more accurate solution 
computed on children cells belonging to level £ + 1. Thus, 
accurate integration of the source term is not important 
for these cells and could even be skipped. From these 
considerations, one may as well evaluate the source term- 
related time step on leaf cells only, where the accuracy 
and stability of the computed solution is essential. This 
trick should speed up the computations by a factor of 
approximately (2)^, thus allowing to take full advantage 
of the refinement offered by the AMR algorithm without 
the time step restriction. Besides, this should not alter 
nor degrade the solution computed during this single hi- 
erarchical integration step as long as the projection step 
precedes the regrid process. 

The proposed modification is expected to be particu- 
larly efficient in those problems where radiative losses are 
stronger in proximity of steep gradients. 

4.3. Parallelization and load balancing 

Both PLUTO and PLUTO-CHOMBO support calcu- 
lations in parallel computing environments through the 
Message Passing Interface (MPI). Since the time evolu- 
tion of the AMR hierarchy is performed recursively, from 
lower to upper levels, each level of refinement is paral- 
lelized independently by distributing its boxes to the set 
of processors. The computation on a single box has no 
internal parallelization. Boxes are assigned to processors 
by balancing the computational load on each processor. 
Currently, the workload of a single box is estimated by 
the number of grid points of the box, considering that 
the integration requires approximately the same amount 
of flops per grid point. This is not strictly true in some 
specific case, e.g in the presence of optically thin radia- 
tive losses, and a strategy to improve the load balance in 
such situations is currently under development. On the 
basis of the box workloads, CHOMBO's load balancer 
uses the Kernighan-Lin algorithm for solving knapsack 
problems. 

CHOMBO weak scaling performance has been thor- 
oughly benchmarked: defining the initial setup by spa- 
tially replicating a 3D hydrodynamical problem propor- 
tionally to the number of CPUs employed, the execution 
time stays constan t with excellent approximation (see 
Van Straalen||2009 for more details). 

To test the parallel performance of PLUTO-CHOMBO 
in real applications, we performed a number of strong 
scaling tests by computing the execution time as a func- 
tion of the number of processors for a given setup. While 
this is a usual benchmark for static grid calculations, in 
the case of AMR computations this diagnostic is strongly 
problem-dependent and difficult to interpret. 

In order to find some practical rule to improve the 
scaling of AMR calculations, we investigated the depen- 
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Figure 6. Density profiles for tiie ID sliock tube problem at 
t = 0.4 with a = TV. Solid lines correspond to admissible ana- 
lytic solutions with (black) and without (red) a compound wave. 
The numerical solution (symbols) is obtained with 400 grid points. 



dency of the parallel performance on some parameters 
characterizing the adaptive grid structure: the maximum 
box size allowed and the number of levels of refinement 
employed using different refinement ratios. As a gen- 
eral rule, the parallel performance deteriorates when the 
number of blocks per level becomes comparable to the 
number of processors or, alternatively, when the ideal 
workload per processor (i.e. the number of grid cells of a 
level divided by the number of CPUs) becomes compara- 
ble to the maximum box size. As we will show, decreasing 
the maximum possible block size can sensibly increase 
the number of boxes of a level and therefore improves 
the parallel performance. On the other hand, using less 
refinement levels with larger refinement ratios to achieve 
the same maximum resolution can lower the execution 
time, reducing the parallel communication volume and 
avoiding the integration of intermediate levels. 

In Section [5] and [6] we will present several parallel scal- 
ing tests and their dependence on the aforementioned 
grid parameters. 

5. MHD TESTS 

In this section we consider a suite of test prob- 
lems specifically designed to assess the performance of 
PLUTO-CHOMBO for classical MHD flows. The selec- 
tion includes one, two and three-dimensional standard 
numerical benchmarks already presented elsewhere in the 
literature as well as applications of astrophysical rele- 
vance. The single-patch integrator adopts the charac- 
teristic tracing step described in section |3.2| with either 
PPM, WENO or linear interpolations carried out in char- 
acteristic variables. 

5.1. Shock tube problems 

The shock tube test problem is a common benchmark 
for an accurate description of both continuous and dis- 
continuous flow features. In the following we consider one 
and three dimensio nal configuration s of s tandard shock 
tubes proposed by Torrilhon ( 2003 ) and Ryu & Jones 
(p95|. 



5.1.1. One- Dimensional Shock tube 



Following | Torrilhon| (12003) we address the capability 
of the AMR scheme to nandle and refine discontinuous 
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Figure 7. Density profiles for tiie ID sliock tube problem at 
t = 0.4 at the vicinity of the compound wave locus, with o = 3. 
The AMR levels vary from to 10 as reported in the legend, explor- 
ing cases of equivalent resolution of 512 (solid), 1024 (dot), 2048 
(dash), 4096 (dot-dash), 8192 (3 dot-dash), 16, 384 (long dash) and 
1,048,576 (solid red) points. At high resolution the solution con- 
verges to the regular one. This figure is analogous to fig. 3 of 
|Promang et al.| ( |2006t . 



features as well as to correctly resolve the non unique- 
ness issu e of MHD Riemann problems in finite vol ume 
schemes ( |Torrilhon|r2003[ [Torrilhon fc Balsara|2004| and 
Lett and right states are given by 



references therein 



(1,0,0,0,1,1,0,1)' 



for xi < , 



I Vfl, = (0.2,0,0,0,l,cos(a),sin(a),G.2)' for xi > , 

(54) 

where V = {p,Vx,Vy,Vz, B^, By, Bz,p) is the ve ctor of 
prirnitive v ariables and F = 2. As discussed by |Torril-| 
hon (2003), for a wide range of initial conditions MHD 
Riemann problems have unique solutions, consisting of 
Lax shocks, contact discontinuities or usual rarefaction 
waves. Nevertheless, there exist certain sets of initial val- 
ues that can result in non unique solutions. When this 
occurs, along with the regular solution arises one that al- 
lows for irregular MHD waves, for example a compound 
wave. A special case where the latter appears is when 
the initial transverse velocity components are zero and 
the magnetic field vectors are a nti-parallel. Such a case 
was noted by Brio & Wu ( 1988 1 and can be reproduced 



by simply choosing a = tt m our initial condition. 

A one dimensional non unique solution is calculated 
using a static grid with 400 zones, x S [—1,1.5]. Left 
and right boundaries are set to outflow and the evolu- 
tion stops at time t = 0.4, before the fast waves reach 
the borders. The resulting density profile is shown in Fig. 
([6|. The solid lines denote the two admissible exact so- 
lutions: the regular (red) and the one containing a com- 
pound wave (black), the latter situated at x ^ —0.24. 
It is clear that the solution obtained with the Godunov- 
type code is the one with the compound wave (symbols) . 

The crucial problematic of this t est occurs when a is 
close to but not exactly equal to tt. Torrilhon (2003) has 
proven that regardless of scheme, the numerical solution 
will erroneously tend to converge to an "irregular" one 
similar to a = tt (pseudo-convergence), even if the initial 
conditions should have a unique, regular solution. This 
pathology can be cured eit her with high order schemes 
(Torrilhon & Balsara 2004) or with a dramatic increase 
in resolution on the region of interest, proving AMR to 



Table 2 

CPU running time for the one-dimensional MHD shock-tube 
using both static and AMR computations. 



Static Run 




AMR Run 




Gain 




Time (s) 


Level 


Ref ratio 


Time (s) 


512 


0.5 





2 


0.5 


1 


1024 


1.9 


1 


2 


1.1 


1.7 


2048 


7.5 


2 


2 


2.3 


3.3 


4096 


31.6 


3 


2 


4.5 


7.0 


8192 


138.0 


4 


2 


8.7 


15.9 


16384 


546.1 


5 


2 


16.9 


32.3 


1048576 


2.237 ■lO'^t 


10 


2 (4) 


1131.1 


1977.8 



Note. — The first and second columns give the number of 
points Nx and corresponding CPU for the static grid run (no 
AMR). The third, fourth and fifth columns give, respectively, 
the number of levels, the refinement ratio and CPU time for the 
AMR run at the equivalent resolution. The last row refers to the 
solid red line of Fig. [t] where a jump ratio of four was introduced 
between levels 6 aiid 7 to reach an equivalent of ~ 10® grid 
points. The last column shows the corresponding gain factor 
calculated as the ratio between static and AMR execution time, 
t CPU time has been inferred from ideal scaling. 

be quite a useful tool. To demonstrate this we choose 
a = 3 for the field's twist. 

Starting from a coarse grid of 512 computational zones, 
we vary the number of refinement levels, with a con- 
secutive jump ratio of two. The 10 level run (solid red 
line) incorporates also a single jump ratio of four be- 
tween the sixth and seventh refinement levels, reaching a 
maximum equivalent resolution of 1,048,576 zones (see 
Fig. [t]). The refinement criterion is s et u pon the variable 
a ^ {Bl + Bl + Bl)/p using Eq. ^ with a thresh- 
old Xr = 0.03, whereas integration is performed using 
PPM with a Roe Riemann solver and a Courant number 
Ca = 0.9. As resolution increases the compound wave 
disentangles a nd the solution converges to the expe cted 



regular form (Torrilhon 2003) Fromang et al 
Table [2] we compare the CPU 



2006 ). In 
running time of the AMR 



runs versus static uniform grid computations at the same 
effective resolution. With 5 and 10 levels of refinement 
(effective resolutions 16,384 and 1,048,576 zones, respec- 
tively) the AMR approach is ^ 32 and ~ 1978 times 
faster than the uniform mesh computation, respectively. 

5.1.2. Three- Dimensional Shock tube 

The second Riemann problem was propo sed by | Ryu & [ 
'Jones' (1995) and later considered also by Toth (2000|; 
Balsara & Shu (.2000) ,, Ga^^ fc Stone (2008) , ^lignone] 



fc 'izeferacosT ( pOTUl ) 
discontinuity is descri 



Mignone et al. An initial 



DCd m terms ot primitive variables 



Vr 



1.08,1.2,0.01,0.5, 



Yr = 1,0,0,0, 



4^^,0.95 
V47r v47r 

for xi <0 

T 

1 ^ 



for xi > 

(55) 

where V = {p,Vi,V2,V3, Bi, B2, B^^p) is the vector of 
primitive variables. The subscript "1" gives the direc- 
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Figure 8. Primitive variable profiles for the 3D shock tube problem at 4 = 0.02 cos a cos 7, along the x direction. Density and thermal 
pressure are plotted in the top panels while vector field components normal ("1") and transverse ("2" and "3") to the initial surface of 
discontinuity are shown in the middle and bottom panels. We show a smaller portion of the domain, x S [—0.25,0.55], in order to emphasize 
the change of resolution by symbol density. 




Figure 9. Closeup of the top-hat feature in the density profile for 
the 3D shock tube problem, along with AMR level structure and 
the mesh. Different colors are used to distinguish grid levels. 



tion perpendicular to the initial surface of discontinuity 
whereas "2" and "3" correspond to the transverse direc- 
tions. We first obtain a one-dimensional solution on the 
domain x € [—0.75, 0.75] using 6144 grid points, stopping 
the computations at t = 0.2. 

In order to test the ability of the AMR scheme to 
maintain the translational invariance and properly re- 



fine the flow discontinuities, the shock tube is rotated in 
a three dimensional Cartesian domain. The coarse level 
of the computational domain consists of [384 x 4 x 4] 
zones and spans [—0.75,0.75] in the x direction while 
y,z £ [0,0.015625]. The rotation angles, 7 around the 
y axis and a around the z axis, are chosen so that 
the planar symmetry is satisfied by an integer shift of 
cells (ria;, n„, ?T-z). Th e rotation matrix can be found in 
Mignone et al. (2010). By choosing tana = — ri/r2 and 
tan [j = tan 7/ cos a = ri/rs, one can show (Gardiner & 



Stone||2008 ) that the three integer shifts 
obey 



I 

r2 



0, 



Uy , Uz must 



(56) 



where cubed cells have been assumed and {ri,r2,r^) = 
(1,2,4) will be used. Computations stop at t = 
0.2 cos a cos 7, once again before the fast waves reach the 
boundaries. We employ 4 refinement levels with con- 
secutive jumps of two, corresponding to an equivalent 
resolution of 6144 x 64 x 64 zones. The refinement cri- 
terion is based on the normalized second derivative of 
X = {\Bx\ + \By\)p with a threshold value Xr = 0.1. Inte- 
gration is done with PPM reconstruction, a Roe Riemann 
solver and a Courant number of C'a = 0.4. 
The primitive variable profiles (symbols) are displayed 
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in Fig. [s] along the x directiorj^ together with the one 
dimensional reference solution in the x € [—0.25,0.55] 
region . In agreement with the solution of|Ryu fc Jones 



(1995), the wave pattern produced consists of a contact 



discontinuity that separates two fast shocks, two slow 
shocks and a pair of rotational discontinuities. A three 
dimensional closeup of the top-hat feature in the density 
profile is shown in Fig. [9] along with AMR levels and 
mesh. The discontinuities are captured correctly, and 
the AMR grid structure respects the plan e symmetry. 
Our results favorably compare with those o f [Gardinerfc 



Stone (20081 ;'Mig none &: Tzeferacos| pOTo| ; [Mignoneet 
aL (2010) and previous similar 2D contigurations. 

The AMR computation took approximately 3 hours 
and 53 minutes on two 2.26 GHz Quad-core Intel Xeon 
processors (8 cores in total). For the sake of comparison, 
we repeated the same computation on a uniform mesh of 
768 X 8 X 8 zones (1/8 of the effective resolution) with 
the static grid version of PLUTO employing « 79 sec- 
onds. Thus, extrapolating from ideal scaling, the com- 
putational cost of the fixed grid calculation is expected 
to increase by a factor 2^^ giving an overall gain of the 
AMR over the uniform grid approach of ^ 23. 

5.2. Advection of a magnetic field loop 

The next test problem considers the two dimensional 
ad vection of a magnet i c field loop. This test, proposed 
by Gardiner & Stone ( 2005 1 , aims to benchmark the 
scheme's dissipative properties and the correct discretiza- 
tion balance of multi-dimensional terms through moni- 
toring the preservation of the initial circular shape of the 
loop. 

As in Gardiner & Stone (2005) and Fromang et al. 
(|2006|), we define the computational domain by x G 



[— 1, IJ and y € [—0.5, 0.5] discretized on a coarse grid of 
64 X 32 grid cells. In the initial condition, both density 
and pressure are uniform and equal to 1, while the veloc- 
ity of the flow is given by v = Vq cos ae^ + Vq sin aey with 
sin a = and cos a = 2/V5. The magnetic 
field is then defined through its magnetic vector potential 
as 

f Ao(R.-r) if r<R, 
if r > R, 



where Aq = IQ-^, R = 0.3, and r = ^/x'^+y^. The 
simulation evolves until t — 2, when the loop has thus 
performed two crossing through the periodic boundaries. 
The test is repeated with 2, 3 and 4 levels of refinement 
(jump ratio of two), resulting to equivalent resolutions 
of [256 X 128], [512 x 256] and [1024 x 512] respectively. 
Refinement is triggered whenever the second deriv ativ e 
error norm of {B^ ^2) X 10^ computed via Eq. (p|, 
exceeds the threshold Xr = 0.1. The integration is car- 
ried out utilizing WENO reconstruction and the Roe Rie- 
mann solver, with Ca = 0.4. 

The temp oral evolution of magnetic energy density is 
seen in Fig. ( 10 ), with 4 levels of refinement. As the field 



loop is transported inside the computational domain, the 
grid structure changes to follow the evolution and retain 



^ Note t hat similar plots were produced in MT and |Mignone| 
|et al.| | |2010| but erroneously labeled along the "rotated direction" 
rattier tnan the x axis. 



the initial circular form. An efficient way to quantita- 
tively measure the diffusive properties of the scheme is 
to monitor the dissip atio n of the magnetic energy. In 
the top panel of Fig. ( 11 ) we plot the normalized mean 
magnetic energy as a function of time. By increasing the 
levels of refinement the dissipation of magnetic energy 
decreases, with < > ranging from ~ 94% to ^ 98% 
of the initial value. In order to quantify the computa- 
tional gain of the AMR scheme we repeat the simula- 
tions with a uniform grid resolved onto as many points 
as the equivalent resolution, without changing the em- 
ployed numerical method. The speed-up is reported in 
the bottom panel of Fig. [TT| 

5.3. Resistive Reconnection 

Magnetic reconnection refers to the process of break- 
ing and reconnection of magnetic field lines with opposite 
directions, accompanied with a conversion of magnetic 
energy into kinetic and thermal energy of the plasma. 
This is believed to be the basic mechanism behind en- 
ergy release during solar flares. The fi rst solution t o the 
problem was given independently by Sweet (1958) and 
Parker (1957), treating it as a two dimensional bound- 
ary layer problem in the laminar limit. 

According to the Sweet-Parker model, the magnetic 
field's convective inflow is balanced by Ohmic diffusion. 
Along with the assumption of continuity, this yields a 
relation between reconnection and plasma parameters. If 
L and 6 are the boundary layer's half length and width 
respectively, we can write the reconnection rate £ as 



£ EE 



"out 



L 



1 



S 



(58) 



With Win and Uout we denote the inflow and outflow 
speeds, into and out of the boundary layer, respectively. 
The Lundquist number for the boundary layer is deflned 
as S = uaL/t], with ua being the Alfven velocity directly 
upstream of the layer, rj the magnetic resistivity and L 
the layer's half length. This dependency of the recon- 
nection rate with the square root of magnetic resistivity 
is called the Sweet -Parker scaling and has been verified 
both numerically ( Biskamg 1986[ [Uzden sky &: Kulsrud 
2000 1 and experimentally 



irtr.||i9y8|) . 



Fmlowing the guidel ines of the GEM Magnetic Re 



2001), the computa- 
Cartesian box, with 



connection Challenge (Birn et al. 
tional domain is a two dimensiona 
X £ [—Lx/2, Lx/2] and y G [—Ly/2,Ly/2] where we 
choose Lr^ = 25.6 and Ly — 12.8. The initial con- 
dition consists of a Harris current sheet: the magnetic 
field configuration is described by B^iy) = i?o tanh(j//A) 
whereas the flow's density is p = posech {y/X) + Poo, 
where A = 0.5, po = 1 and poo = 0.2. The flow's thermal 
pressure is deduced assuming equilibrium with magnetic 
pressure, P = Bq/2 = 0.5. The initial magnetic field 
components arc perturbed via 



dB^ = -'^o{^/Ly)sm{Try/Ly)cos{2Trx/L^) 
dBy — 'i'o{2'K/Lx)sin{2'Kx/Lx) cos{ny /Ly). 



(59) 



where ^l^o = 0.1. The coarse grid consists of [64 x 32] 
points and additional levels of refinement are triggered 
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Figure 10. Magnetic energy density for the 2D field loop problem at f = 0, 0.7, 1.4, 2. Overplotted are the refinement levels. 
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Figure 11. Upper panel: Normalized magnetic energy for the 
field loop advection problem as a function of time. Lower panel: 
computational time as a function of equivalent resolution. 



using the following criterion based on the current density: 



AyBj. 



> 



Xn 



\A,By\ + \AyB,\+^/C I +{1-0 



Xn 



(60) 



where A^rBy and AyB. 



are the undivided central differ- 
ences of By and B^ in the x and y direction, respectively, 
and ^ — Ax'^/Ax^ > 1 is the ratio of grid spacings be- 
tween the base (0) and current level (l). The threshold 
values Xmin — 0.2 and Xmax = 0.37 are chosen in such 
a way that refinement becomes increasingly harder for 
higher levels. We perform test cases using either 5 levels 
of refinement with a consecutive jump ratio of two or 3 
levels with jump ratio 2:4:4 reaching, in both cases, an 
equivalent resolution of 2048 x 1024 mesh points. Bound- 
aries are periodic in the x direction, whereas perfectly 
conducting boundary walls are set at y = zLLy/2. We 
follow the computations until t — 100 using PPM re- 
construction with a Roe Ri ema nn solver and a Courant 
number Cq = 0.8. In Fig. 12 we display the temporal 
evolution of the current density for a case where the uni- 
form resistivity is set to 77 = 8 • 10~^. A reconnection 
layer is created in the center o f the domain, which pre- 
disposes resistiverecomiectiOT In agree- 
ment with Reynolds et al. ( 2005 1 the maximum value of 
the current density decreases with time, (fig. 4 of that 
study). As seen in the first two panels (t = 0, 50), the re- 
finement criterion is adequate to capture correctly both 
the boundary layer and the borders of the magnetic is- 
land structure. In the rightmost panel (t = 100) we also 
draw sample magnetic field lines to better visualize the 
reconnection region. 

In order to compare our numerical results with the- 
ory, we repeated the computation varying the value of 
the magnetic resistivity ry. For small values of resistivity 
(large Lundquist numbers S), the boundary lay er is elon- 
gated and presents large aspect ratios A = L/6. Biskamp 



( |1986| reports that for A beyond the critical value ot 
Acrit — 100 the boundary layer is tearing unstable, lim- 
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Figure 12. Upper row: temporal evolution of the current density and magnetic field lines for the resistive reconnection problem, with 
ri = 810-3. Snapshots refer to t = 0, 50, 100. Lower row: Pressure profiles for various values of resistivity rj, along with the AMR level 
structure at t = 50. The refinement strategy consists of 3 levels with a jump ratio of 2 : 4 : 4 (equivalent resolution of 2048 X 1024 mesh 
points). 
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Figure 13. Time average of <5/L, analogous to the magnetic recon- 
nection rate £, as a, function of resistivity. Symbols represent the 
actual data, whereas over-plotted (dashed line) is the Sweet-Parker 
scaling ~ y'r;, along with a best fit (solid line). The numerical re- 
sults arc in agreement with the theoretical scaling. 
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iting the resistivity range in which the Sweet-Parker re- 
connection model operates. Since in the Sun's corona 
S can reach values of ^ 10^"* >> S'max, sec ondary is 



land format ion must be taken into account ( Cassak & 
Drake 2009). In this context, ensuring that we respect 
Biskamp's stability criterion, we calculate the temporal 
average oi 5/L, analogous to the reconnection rate £, for 
various rj and reproduce the Sweet-Parker scaling (Fig. 
13). The boundary layer's half -width (S) and -length 
(IT) are estimated from the e-folding distance of the peak 
of the electric current, while the AMR scheme allows us 
to economically resolve the layer's thickness with enough 
grid points. 

Parallel performance for this problem is shown in Fig 



14 



where we plot the speedup 5* = Ti/Tmc-pu as a 
Tunction of the number of processors A^cpu for the 3- 
and 5-level AMR computations as well as for fixed uni- 
form grid runs carried out at the equivalent resolution of 



Figure 14. Parallel speedup at t = 50 as a function of the number 
of processors (A'cpu) for th^ resistive reconnection problem. The 
different lines refer to the execution times obtained with 5 levels 
of refinement with consecutive jumps of 2 (red squares), 3 levels 
with jump ratios 2:4:4 (green crosses) and a fixed uniform grid 
with 2048 X 1024 zones (black plus signs). The dotted line gives 
the ideal scaling whereas the number of blocks on the finest level 
at the end of integration is reported above each curve. 

2048 X 1024 zones. Here Ti is the same reference con- 
stant for all calculations and equal to the (inferred) run- 
ning time of the single processor static mesh computation 
while TjVcpu is the execution time measured with A'^cpu 
processors. The scaling reveals an efficiency (defined as 
S/Nqpu) larger than 0.8 for less than 256 processors with 
the 3- and 5-level computations being, respectively, 8 to 
9 and 4 to 5 times faster than the fixed grid approach. 
The number of blocks on the finest level is maximum 
at the end of integration and is slightly larger for the 
3-level run (1058 vs. 835). This result indicates that us- 
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ing fewer levels of refinement with larger grid ratios can 
be more efficient than introducing more levels with con- 
secutive jumps of 2, most likely because of the reduced 
integration cost due to the missing intermediate levels 
and the decreased overhead associated with coarse-fine 
level communication and grid generation process. Effi- 
ciency quickly drops when the number of CPU tends to 
become, within a factor between 2 and 3, comparable to 
the number of blocks. 

5.4. Current Sheet 



The curren t sheet problem proposed by [Gardiner fc 
Stone (2005) and later considered by Fromang et al. 
I pOOS I in the AMR context is particularly sensitive to nu- 
merical diffusion. The test problem follows the evolution 
of two current sheets, initialized through a discontinuous 
magnetic field configuration. Driven solely by numerical 
resistivity, reconnection processes take place, making the 
resulting solution highly susceptible to grid resolution. 

The initial condition is discretized onto a Cartesian two 
dimensional grid x, y € [0, 2], with 64 x 64 zones at the 
coarse level. The fluid has uniform density p — 1 and 
thermal pressure P = 0.1. Its bulk flow velocity v is set 
to zero, allowing only for a small perturbation in Vx = 
Vq sin(7r?/), where vq = 0.1. The initial magnetic field 
has only one non-vanishing component in the vertical 
direction. 



By - 



-Bo if |a;- 1| < 0.5, 
Bq otherwise , 



(61) 



where Bq = I, resulting in a magnetically dominated 
configuration. Boundaries are periodic and the integra- 
tion terminates at t = 4. We activate refinement when- 
ever th e m aximum between the two error norms (given 



by Eq. 50 1 computed with the specific internal energy 



and the y component of the magnetic field exceeds the 
threshold value = 0.2. In order to filter out noise 
between magnetic island we set e = 0.05. We allow for 
4 levels of refinement and carry out integration with the 
Roe Riemann solver, the improved WENO reconstruc- 
tion and a CFL number of 0.6. 

In Fig. [15] we show temporal snapshots of pressure 
profiles for t = 0.5, 1.0, 1.5, 2, 3.5 and 4.0, along with 
sample magnetic field lines. Since no resistivity has 
been specified, the elongated current sheets are prone 
to tearing instability as numerical resistivity is small 
with respect to Biskamp's criterion (jBiskamp 1986). Sec- 
ondary islands, "plasmoids", promptly form and prop- 
agate parallel to the field in the y direction. As re- 
connection occurs, the field is dissipated and its mag- 
netic energy is transformed into thermal energy, driving 
Alfven and compressional waves which further seed re - 
connection ( Gardiner fc Stone|2005 Lee fc Deane|2009" ). 
Due to the dependency of numerical resistivity on ttie 
field topology, reconnection events are most probable at 
the nodal points of Vx- At the later stages of evolu- 
tion, the plasmoids eventually merge in proximity of the 
anti-nodes of the transverse speed, forming four l arger 
islands, in agreemen t with the results presented in |Gar-| 
diner & Stone ( 2005 1. A closeup on the bottom left island 
is shown in l''igTjl6[ at the end of integration. The refine- 
ment criterion on the second derivative of th ermal pres- 
sure, as suggested by Fromang et al. (2006), efflciently 



captures the island features of the solution. 

5.5. Three-dimensional Rayleigh Taylor Instability 

In this example we consider the dynamical evolution 
of two fluids of different densities initially in hydro- 
static equilibrium. If the lighter fluid is supporting 
the heavier fluid against gravity, the configuration is 
known to be subject to the Rayleigh- Taylor instability 
(RTI). Our computational domain is the box spanned by 
x,z e [-1/2, 1/2], y G [-3/2, 1/2] with gravity pointing 
in the negative y direction, i.e. g = (0,-1,0). The two 
fluids are separated by an interface initially lying in the 
xz plane at y = 0, with the heavy fluid p = ph at y > 
and p — pi at y < 0. In the current example we em- 
ploy ph — 4:, pi — I and specify the pressure from the 
hydrostatic equilibrium condition: 



piy) 



100 

7 



py, 



(62) 



so that the sound crossing time in the light fluid is 0.1. 
The seed for instability is set by perturbing the velocity 
fleld at the center of the interface: 



exp [— (5r)^] 
10cosh(10y)2 



(63) 



where r — \/x^ + z'^. We assume a constant magnetic 
fleld B = {Bx,0,0) parallel to the interface and ori- 
ented in the x direction with different field strengths, 
Bx = (hydro limit), Bx — 0.2Bc (moderate field) and 
Bx = 0.6Bc (strong field). Here B^ = \/{ph - Pi)gL 
is the critical field value above which instabilities are 



suppressed ( Stone k Gardiner||2007 1. We use the PPM 
method witn the Roe Riemann solver and a Courant 
number Ca — 0.45. The base grid has 16 x 32 x 16 
cells and we perform two sets of computations using i) 
4 refinement levels with a grid spacing ratio of 2 and ii) 
2 refinement levels with a grid jump of 4, in both cases 
achieving the same effective resolution (256 x 512 x 256). 
Refinement is triggered using the second derivative error 
norm of density and a threshold value Xrof = 0.5. Pe- 
riodic boundary conditions are imposed at the x and z 
boundaries while fixed boundaries are set at y = —3/2 
and y = 1/2. 

Results for the un-magnetized, weakly and stro ngly 
magnetized cases are shown at < = 3 in Fig. [TtI 
In all cases, we observe the development of a central 
mushroom-shaped finger. Secondary instabilities due to 
Kelvin-Helmholtz modes develop on the side and small- 
scale structures are gradually suppressed in the direc- 
tion of the field as the strength is i ncreas ed. Indeed, as 



pointed out by Stone &: Gardiner] (]2007|), the presence 
of a uniform magnetic field has the efl'ects of reducing 
the growth rate of modes parallel to it although the in- 
terface still remains Rayleigh- Taylor unstable in the per- 
pendicular direction due to interchange modes. As a net 
effect, the evolution becomes increasingly anisotropic as 
the field strengthens. 

Parallel performance is plotted, for the moderate field 
case, in Fig. [18] where we measure the speedup factors 
of the 4- and 2- level computations versus the number 
of CPUs. Here the speedup is defined by S* = Ti/Tncpu 
where Ti is the inferred running time relative to the 4- 
level computation on a single-processor. The two-level 



AMR Computations with the PLUTO Code 



17 




0.5 1.0 1.5 

X 



0.5 1.0 1.5 

X 



Figure 15. Time evolution of tlie pressure profiles along with sample magnetic field lines, for the current sheet problem. Temporal 
snapshots refer to t = 0.5, 1, 1.5, 2 (upper four) and t = 2.5, 3, 3.5, 4 (lower four). 



pressure 




action with the interstellar medium. Energy, mass and 
momentum exchange leading to cloud-crushing strongly 
depends on the orientation of the magnetic field and the 
resulting an isotropy of thermal co nduction. 

Following Orlando et al. (2008), we consider the two- 
dimensional ClartesIan~3omam X G [—4,4], y e [—1.4 x 
6.6] with the shock front propagating in the positive y 
direction and initially located at y = —1. Ahead of the 
shock, for y > — 1, the hydrogen number density ni has 
a radial distribution of the form: 



ni 



cosh [cr{r/rc\)° 



(64) 



^ is the hydrogen number density at 
Ua = O.lric is the ambient number 



Figure 16. Closeup of the bottom left island at t = 4 for the cur- 
rent sheet problem. Pressure contours, along with the refinement 
levels and grid are shown. 

case shows improved performance over the four-level cal- 
culation both in terms of CPU cost as well as parallel effi- 
ciency S/Nqptj . Indeed the relative gain between the two 
cases approaches a factor of 2 for an increasing number 
of CPUs. Likewise, the efhciency of the fewer level case 
remains above > 0.8 up to 512 processors and drops to 
~ 0.61 (versus ~ 0.44 for the 4-level case) at the largest 
number of employed cores (1024) which is less than the 
number of blocks on the finest grid. 

5.6. Two-dimensional Shock-Cloud Interaction 

The shock-cloud interaction problem has been exten- 
sively studied and used as a standard benchmark for the 
vali dation of MHD schemes and inter-code comparisons 



where Uc — 1 cni^ 
the cloud's center 
density, — Ipc is the cloud's radius {a = 10) and r 
is the radial distance from the center of the cloud. The 
ambient medium has a uniform temperature Ta = 10^ K 
in pressure equilibrium with the cloud, with a thermal 
pressure of pi = 2kBnaTa (we assume a fully ionized 
gas). Downstream of the shock, density and transverse 
components of the magnetic field are compressed by a 
factor of n-ijnx — (F H- l)/(r — 1) while pressure and 
normal velocity are given by 



I 2 2kBTs 
F — 1 firuH 



(65) 



seepai fc Woodward|1994[ |Balsara |2001[|Lee k Deane 



2009 and reference therein). In an astrophysical con 



where /i = 1.26 is the inverse of the hydrogen mass frac- 
tion and Ts — 4.7 • 10^ K is the post-shock temperature. 
The normal component of the magnetic field (By) re- 
mains continuous through the shock front. 

We perform two sets of simulations with a magnetic 
field strength of |B| — 1.31 fj,G, initially parallel (case 
Bx4) or perpendicular (case By 4) to the shock front . 



text, it adresses the fundamental issue of the complex 
morphology of supernova remnants as well as their inter- 



Adopting the same notations as Orlando et al 
we solve in each case the MHD equations with ('. 



(2008) 
t^NJ^d 



without (NN) thermal conduction effects for a total of 4 
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Figure 17. Three-dimensional view of the Rayleigh- Taylor instability problem showing density at t = 3 for the un-magnetized (left panel), 
weakly magnetized (middle panel) and strongly magnetized case (right panel) using 4 levels of refinement. In each panel we also show a 
sliced boxed emphasizing the grid refinement ratios. 
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Figure 18. Parallel speedup versus number of processors A^qpu 
for the 3D Rayleigh- Taylor problem computed with 2 (green 
crosses) and 4 (red squares) refinement levels between t = 2 and 
t = 2.25. The dotted line gives the ideal scaling. The number of 
blocks on the finest level at t = 2 is printed above and below the 
corresponding curves. 

cases. Thermal conductivity coefficients along and across 
the magnetic field lines (see Eq. |6]) are given, in c.g.s 
units, by 



9.2-10"^T^/^ fci =5.4-10"^^ 



H 
TB2 



(66) 



The resolution of the base grid is set to 32'^ points and 
5 levels of refinement are employed, yielding an equiv- 



tracing (Section 3.2) together with PPM reconstruction 
(Eq. 33 ) and the HLLD Riemann solver are chosen to ad- 
vance the equations in time. Open boundary conditions 
are applied on all sides, except at the lower y boundary 
where we keep constant inflow values. The CFL number 
is Ca = 0.8. Refinement is triggered upon density with 
a threshold Xrcf = 0.15. 

In Fig. [T9| and [20| we show the cloud evolution for the 
four different cases at three different times t = 1,2,3 
(in units of the cloud crushing time Tcc = 5.4 • 10'^ yrs). 
Our results a re in excellent agreement with those of |0r-| 
lando et al. (20081 confirming the efficiency of thermal 
conduction m suppressing the development of hydrody- 
namic instabilities at the cloud borders. The topology 
of the magnetic field can be quite effective in reducing 
the efficiency of thermal conduction. For a field parallel 
to the shock front (TN-Bx4), the cloud's expansion and 
evaporation are strongly limited by the confining effect 
of the enveloping magnetic field. In the case of a per- 
pendicular field (TN-By4), on the contrary, thermal ex- 
change between the cloud and its surroundings becomes 
more efficient in the upwind direction of the cloud and 
promotes the gradual heating of the core and the conse- 
quent evaporation in a few dynamical timescales. In this 
case, heat conduction is strongly suppressed laterally by 
the presence of a predominantly vertical magnetic field. 

5.7. Radiative Shocks in Stellar Jets 

Radiative jet models with a variable ejection velocity 
are commonly adopted to reproduce the knotty struc- 
ture observed in coUimated outflows from young stel- 
lar objects. The time variability leads to the formation 
of a chain of perturbations that eventually steepen into 
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Figure 19. Top panel: density maps (gr/cm^ in Log scale) for the shock-cloud interaction at t = 1,2,3 for a magnetic field initially 
parallel to the shock interface (Bx4). Here the time unit is defined as the cloud crushing time and equal to 5,400yrs. Bottom panel: 
magnetic field strength (in units of 10~® G) with the distribution of patches at levels 3 and 4 superimposed. The left and right half of each 
panel shows the results computed with the PPM method without (NN) and with (TN) thermal conduction, respectively. The resolution of 
the base grid is 32^ and 5 levels of refinement are used. 
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Figure 21. Radiative jet structure at t ft: 64.3 yrs. On the left we display a full view of the temperature distribution (right) and AMR 
block structure (left). The smaller panels to the right show closeup views of number density (left half, in cm~'^), temperature (right half, in 
10* K), neutral hydrogen (left half) and magnetization i3^/(2p) (right half) for the second (top panels) and third (bottom panels) shock, 
respectively. The AMR block structure is overplotted in the right half of each panel. 
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Figure 22. One dimensional enlarged cuts at r 0.1133rj across 
the middle shock in the radiative jet model at t = 64.3 yrs. Tem- 
perature density distributions are plotted in the top panel while the 
first ionization stage of O, H, N and neutral Sulphur are plotted in 
the bottom panel. 



pairs of forward/reverse s hocks (the "internal working 
surfaces" see, for instance de CoUe fc Raga]|2006[ |Raga| 



et al. 12007 Te§ileanu et al. 2009 ) which are neld respon 
3 to 



sible tor the typical spectral signature of these objects. 
While these internal working surfaces travel down the 
jet, the emission from post-shocked regions occurs on a 
much smaller spatial and temporal scale thus posing a 
considerable challenge even for AMR codes. 

As an example application, we solve the MHD equa- 
tions coupled to the chemical network described in 



Te§ileanu et al. (2008), evolving the time-dependent ion- 
ization and collisionally-excited radiative losses of the 
most important atomic/ionic species: H, He, C, N, O, 
Ne and S. While hydrogen and helium can be at most 
singly-ionized, we only consider the first three (instead 
of five) ionization stages of the other elements which suf- 
fices for the range of temperatures and densities con- 
sidered here. This amounts to a total of 19 additional 
non-homogeneous continuity/rate equations such as Eq. 

& 

I'he initial condition consists of an axisymmetric su- 
personic coUimated beam in cylindrical coordinates (r, z) 
in equilibrium with a stationary ambient medium. Den- 
sity and axial velocity can be freely prescribed as 



p{r) 



Pj - Pa 
cosh(r*^/r^) 



cosh(r''/r^) 



(67) 

where pj and pa — are the jet and ambient mass 

density, vj = 200km/s is the jet velocity and rj = 
2 • 10"'^^ cm is the radius. The jet density is given by 
Pj 



jfijiria where 



— in3 



10 is the total particle num- 



ber density, is the mean molecular weight and nia is 
the atomic mass unit. Both the initial jet cross section 
and the environment are neutral with the exception of C 
and S which are singly ionized. The steady state results 
from a radial balance between the Lorentz and pressure 



forces and has to satisfy the time-independent momen- 
tum equation: 



dp 
dr 



1 1 d{rB^ 



2r2 



dr 



where, for simplicity, we neglect rotations and assume a 
purely azimuthal magnetic field of the form 



B^{r) = ^Vl-cxp [-(r/a)4] 



(69) 



Here a = 0.9rj is the magnetization radius and is a 
constant that depends on the jet and ambient tempera- 
tures. Direct integration of Eq. (68) yields the equilib- 
rium profile 



p{r) = pj 



IBl^eriirya^) 
2 a2 



(70) 



where pj = pjksTj / [pjiria) is the jet pressure on the 
axis, jij is the mean molecular weight, ma is the atomic 
mass unit, ks is the Boltzmann constant and Tj is the 
jet temperatu re. The value of B^. is recovered by solv- 
ing equation (70 1 in the r — > cx) limit given the jet and 
ambient temperatures Tj = 5 • 10"^ K and = 10'^ K, 
respectively. O ur choice of param eters is similar to the 

(20071). 



one adopted by Raga et al. 
The computational domain 



defined by r/rj e 



[0, 12.5], z/rj € [0, 25] with axisymmetric boundary con- 
ditions holding at r = 0. At 2: = we keep the flow 
variables constant and equal to the equilibrium solution 
and add a sinusoidal time-variability of the jet veloc- 
ity with a period of 20 yrs and an amplitude of 40km/s 
around the mean value. Free-outflow is assumed on the 
remaining sides. We integrate the equation s us ing linear 
reconstruction with the MC limiter, Eq. (34), and the 



HLL C Riem ann solver. The time step is computed from 
Eq. (23 24) using a Courant number Ca = 0.6, a rela- 
tive tolerance ec — .02 and following the considerations 
given in Section 4.2 Shock-adaptive hybrid integration 



is used by locally switching to the more diffusive MinMod 
limiter (Eq. 35 ) and the HLL Riemann solver, according 
to the mechanism outlined in Appendix [B| 

The base grid consists of 128 x 256 zones and 5 addi- 
tional levels with grid refinement jumps of 2 : 2 : 2 : 4 : 4 
are used. At the effective resolution of 16384 x 32768 
zones the minimum length scale that can be resolved 
corresponds to 1.5 26 ■ 10^^ cm (s a 0.1 AU), the same as 



model M3 of Raga et al. (2007). Zones are tagged for 
refinement when the second derivative error norm of den- 
sity exceeds the threshold value Xrcf = 0.15. In order to 
enhance the resolution of strongly emitting shocked re- 
gions, we prevent levels higher than 2 to be created if 
the temperature does not exceed 5 • 10'' x £ where £ is the 
level number. 

The results are shown, after « 63.4 years, in Figure 
[2T] where the steepening of the perturbations leads to 
the formation of forward/reverse shock pairs. Radiative 
losses become strongly enhanced behind the shock fronts 
where temperature attains larger values and the com- 
pression is large. This is best illustrated in the closeup 
views of Fig [21] where we display density, temperature, 
hydrogen fraction and magnetization for the second (up- 
per panel) and third shocks (lower panel) , respectively lo- 
cated at z ^ 14.2 and z sa 7.34 (in units of the jet radius). 



22 



Mignone et al. 





a 








6 

4 
2 

0. 


35 0.87 0.90 0. 


32 













n 












0-1 








0.0 
-0.1 
-0,2 








-0-3 








-0.4 
-0.5 


ai 






0,85 0,87 0,90 0, 


32 




0.2 


0.4 0.6 


0.8 



Figure 23. Density (top left), longitudinal velocity (top right), y 
components of magnetic field and velocity (bottom left and right) 
for the one-dimensional relativistic magnetized shock tube prob- 
lem at t=0.4. We employ a base grid of 400 cells with 5 levels of 
refinement v^ith consecutive jump ratios of 4 : 2 : 2 : 2 : 2 yield- 
ing an equivalent resolution of 25600 zones. The grid hierarchy is 
shown in the top left panel (dashed red line) while magnifications 
of the thin shell are plotted, for each quantity, using symbols in 
the interior plot windows. 

Owing to a much shorter cooUng length, tc ~ p/A, a thin 
radiative layer forms the size of which, depending on local 
temperature and density values, becomes much smaller 
than the typical advection scale. An extremely narrow 
one-dimensional cut shows, in Figure |22[ the profiles of 
temperature, density and ionization fractions across the 
central radiative shock, resolved at the largest refinement 
level. Immediately behind the shock wave, a flat tran- 
sition region with constant density and temperature is 
captured on ^ 6 grid points and has a very small thick- 
ness ~ 0.005rj. A radiatively cooled layer follows behind 
where temperature drops and the gas reaches the largest 
ionization degree. Once the gas cools below 10'* K, radia- 
tive losses become negligible and the gas is accumulated 
in a cold dense adiabatic layer, A proper resolution of 
these thin layers is thus crucial for correct and accurate 
predictions of emission lines and int ensity ratios in stel- 
lar jet models (Te§ileanu et al,|2011 ). Such a challenging 
computational problem can be tackled only by means of 
adaptive grid techniques, 

6, RELATIVISTIC MHD TESTS 

In this section we apply PLUTO-CHOMBO to test 
problems involving relativistic magnetized flows in one, 
two and three dimensions. Both standard numeri- 
cal benchmarks and applications will be considered. 
By default, the conservative MUSCL-Hancock predictor 
scheme (Eq. 43 1 together with linear reconstruction on 
primitive variables are used during the computation of 
the normal predictors, 

6,1. One dimensional Shock-Tube 

Our flrst test consists of an initial discontinuity located 
at a; = 0,5 separating two regions of fluids characterized 
by 



(By, B,, p) 



(7, 7, 103) 
(0,7, 0.7, 0,1) 



for 
for 



X < 0,5 
X > 0,5 



Table 3 

CPU performance for the one-dimensional RMHD 
shock-tube. 



Static Run 




AMR Run 




Speedup 


Nx Time (s) 


Level 


Ref ratio 


Time (s) 


1600 6,1 


1 


4 


1.8 


3.4 


3200 29,9 


2 


2 


3.8 


7.9 


6400 124.5 


3 


2 


8.2 


15,2 


12800 502.3 


4 


2 


18.1 


27,8 


25600 2076.1 


5 


2 


41.2 


50.4 



Note. — The first and second columns give the number of 
points Nx and corresponding CPU for the static grid run (no 
AMR). The third, fourth and fifth columns give, respectively, 
the number of levels, the refinement ratio and CPU time for 
the AMR run at the equivalent resolution. The last column 
shows the corresponding speedup factor. 



(see Balsar a |200l| [Mignone fc Bodo|2006 

et al. 2008iaiid reierences therein), I'he t 



van der Hoist 
), 'I'he fluid is initially 
1 and the longitudinal 



iet ai._,zuuo| £ 

at rest with uniform density p 
magnetic fleld is = 10. 

We solve the equations of relativistic MHD with the 
ideal gas law ( F = 5/3) using t he 5- w ave HLLD Rie- 
mann solver of [Mignone et al,[ mom and Ca = 0,6. 
The computational domain [0, 1] is discretized using 5 
levels of refinement with consecutive grid jump ratios 
of 4 : 2 : 2 : 2 : 2 starting from a base grid (level 
0) of 400 grid zones (yielding an effective resolution of 
25600 zones). Refinement is trig gere d upon the vari- 
able a ^ {B^ + Bl)/D using Eq. |50|) with a threshold 
~Xr = 0,1, At t = 0,4 the resulting wave pattern (Fig. 
23) is comprised of two left-going rarefaction fans (fast 
and slow) and two right-going slow and fast shocks. The 
presence of magnetic fields makes the problem particu- 
larly challenging since the contact wave, slow and fast 
shocks propagate very close to each other, resulting in a 
thin relativistic blast shell between x « 0,88 and x « 0,9. 

In Table [3| we compare, for different resolutions, the 
CPU timing obtained with the static grid version of the 
code versus the AMR implementation at the same effec- 
tive resolution, starting from a base grid of 400 zones. 
In the unigrid computations, halving the mesh size im- 
plies approximately a factor of four in the total running 
time whereas only a factor of two in the AMR approach. 
At the large resolution employed here, the overall gain is 
approximately 50. This example confirms that the res- 
olution of complex wave patterns in RMHD can largely 
benefit from the use of adaptively refined grids, 

6,2. Inclined Generic Alfven Test 
The generic Alfven test ( Giacomazzo fc Rezzolla|2006 



Mignone et al,|2d09 | consists ot the tollowmg non planar 
initial discontinuity 

f Yl = (1,0,0,3,0,4,1,6,2,5)^ for 2:i<l/2, 



V, 



(0,9, 0, 0, 0, 1, 5, 2, 5,3)^ for xi > 1/2 , 



(71) 



, , (72) 

where, as in Section [5,1.1[ V = 
(p^vi,V2,V3, Bi, B2, B^jp) is given in the frame of 
reference aligned with the direction of motion xi. The 
ideal equation of state (10) with F = 5/3 is adopted. 
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Figure 24. Horizontal cuts for the rotated inclined Alfven test at t = 0.4/\/2 (symbols) and for the 1-D reference solution at t = 0.4 
(solid line). The top row shows, from left to right, proper density, gas pressure and Lorentz factor. In the middle and bottom rows we plot 
the components of velocity and magnetic field normal {vi and Bi) and transverse {v2, f3 and B2, Bs) to the surface of discontinuity. 



at the top and bottom boundaries where, for any flow 
quantity q, we set q{i,j) = q{i ± 1, j =F 1). 



Sigma 



AMR Level 




Figure 25. A closeup of the central region in the inclined Alfven 
test at t = 0.4/\/2 showing a = B^/D with the overplotted block 
distribution. The base grid has 64 X 4 zones and 5 levels of refine- 
ment are used. 



Here we consider a two-dimensional version by rotat- 
ing the discontinuity front by 7r/4 with respect to the 
mesh and following the evolution imtil t — 0.4cos7r/4. 
The test is run on a coarse grid of 64 x 4 zones cov- 
ering the domain x € [0,1], y € [—1/32,1/32] with 6 
additional levels of refinement resulting in an effective 
resolution of 4096 x 256 zones (a factor of 2 in resolution 
is used between levels). In order to trigger refinement 
across jumps of different natu re, we define the quantity 
<T = 'B^/D together with Eq. (501 and a threshold value 
Xr = 0.03. The Courant nurnBer is C'a = 0.6 and the 
HLLD Riemann solver is used throughout the computa- 
tion. Boundary conditions assume zero-gradients at the 
X boundaries while translational invariance is imposed 
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The breaking of the discontinuity, shown in Fig 
leads to the formation of seven waves including a IctF 
going fast rarefaction, a left-going Alfven discontinu- 
ity, a left-going slow shock, a tangential discontinuity, 
a right-going slow shock, a right-going Alfven discon- 
tinuity and a right-going fast shock. Our results indi- 
cate that refined regions are created across discontinuous 
fronts which are correctly followed and resolved, in ex- 
cellent agreement with the reference solution (obtained 
on a fine one-dimensional grid). Note also that the lon- 
gitudinal component of the m agn etic field, Bi, does not 
present spurious jumps (Fig. |24[ ) and shows very small 
departures from the expected constant value. In the cen- 
tral region, for 0.4 < x < 0.6, rotational discontinuities 
and slow shocks are moving slowly and very adjacent to 
each other thus demanding relatively high resolution to 
be correctly captured. With the given prescription, grid 
adaptation efficiently provides adequate resolution where 
needed. This is best shown in the closeup of the central 
region, Fig. [25) showing a together with the AMR block 
structure. 

A comparison of the performance between static and 
adaptive grid computations is reported in Table l4] using 
different resolutions and mesh refinements. At the reso- 
lution employed here, the gain is approximately a factor 
of 9. 



24 



Mignone et al. 



Table 4 

CPU performance for the two-dimensional inclined generic Alfvcn test 



Static Run 



AMR Run 



Speedup 





Time (s) 


Level 


Ref ratio 


#Blocks 


Time (s) 




128 


1.8 


1 


2 


4 


1.1 


1.6 


256 


10.8 


2 


2 


14 


6.1 


1.8 


512 


75.2 


3 


2 


44 


33.0 


2.3 


1024 


557.8 


4 


2 


126 


185.9 


3.0 


2048 


4390.5 


5 


2 


306 


900.0 


4.9 


4096 


35250.0 


6 


2 


761 


4346.9 


8.1 



6.3. Relativistic Rotor Problem 



Th e two-dimensional rotor problem 



2003| [van der Hoist et al.|[2008} [Dumbser ef al. 2001 



Del Zanna et al 



consists ot a rapidly spinning disk embedded in a uni 
form background medium (p = p = 1) threaded by a 
constant magnetic field B = (f,0, 0). Inside the disk, 
centered at the origin with radius R = O.I, the density is 
p = 10 and the velocity is prescribed as v = uj{—y, x, 0) 
where a; = 9.95 is the angular frequency of rotation. The 
computational domain is the square x^y G 5] with 
outflow (i.e. zero-gradient) boundary conditions. Com- 
putations are performed on a base grid with 64^ grid 
points and 6 levels of refinement using the HLL Riemann 
solver and a CFL number Ca = 0.5. Zones are tagged for 
refinement whenever the sec ond derivative error norm of 
£/{fry), computed with (50), exceeds Xr = 0.15. The ef- 
fective resolution amounts therefore to 4096^ zones, the 
largest employed so far (to the extent of our knowledge) 
for this particular problem. 

Results are shown at t = 0.4 in Fig. [26] where one can 
see an emerging flow structure enclosed by a circular fast 
forward shock (r w 0.425) traveling into the surrounding 
medium. An inward fast shock bounds the innermost 
oval region where density has been depleted to lower val- 
ues. The presence of the magnetic field slows down the 
rotor and the maximum Lorentz factor decreases from 
10 to « 2.2. The numerical solution preserves the ini- 
tial point symmetry and density corrugations, present in 
lower resolution runs, seem drastically reduce d at this 
resolution, i n acc ordance with the findings of |van der| 



Hoist et al. (2008). Divergence errors are mostly gener 
ated in proximity of the outer fast shock and are damped 
while being transported out of the domain at the speed 
of light in the GLM formulation. We monitor the growth 
of such errors by plotting the volume-average of |V • B| 



as a function of time showing (bottom panel in Fig. 26 ) 
that the growth of monopole errors is limited to the same 
values of van der Hoist et al.| f2008 ) (the factor of 4 comes 
from the fact that our computational domain is twice as 
small) . 

Parallel performance for this particular test is shown in 
Fig. [27| p lotting the speedup, defined as S" = ST^/Tat^p^, 
where is the execution time of the uniform grid run 
performed at the same effective resolution (4096^) on 8 
processors while T/v^pu is the running time measured 
with A^cpu processors. AMR scaling has been quanti- 
fied for computations using maximum patch sizes of 16 
and 32 grid points resulting, respectively, in II92I and 
6735 total number of blocks at the end of integration. 



Approximately half of it belongs to the finest level, as 
reported in Fig 27 In general, the AMR approach offers 
a 5 to 10 speedup gain in terms of execution time over the 
tmiform, static grid approach. For iVcpu < 256 the par- 
allel efficiency, measured as S/Nqy>Vi is larger than 0.7 
while it tends to be reduced for more processors. Compu- 
tations carried with fewer blocks (i.e. larger block sizes) 
tend to be more efficient when fewer CPUs are in use 
since inter-processor communications are reduced. How- 
ever, this cost becomes sensibly larger at 512 (or more) 
processors resulting in a loss of efficiency. 



Note. — The base grid for the AMR computation is 64 X 4 zones. 
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6.4. Cylindrical Blast Wave 

Strong symmetric explosions in highly magnetized 
environments can become rather arduous benchmarks 
probing the code's abilit y to evolv e strongly magn etized 
sh ocks (see, for instance iKomissa rov 1999 

'ar'2005'; 'Mignone fc Bodo"!jOO^! |Dei Zanr 

,Beckwith fc Stone 20 1 1 ) . Failures may lead to unphysi 
cal densities or pressures whenever the scheme does not 
introduce adequate dissipation across oblique discontinu- 
ities and/or if the divergence- free condition is not prop- 
erly controlled. 

The two dimensional setup considered here consists 
of the Cartesian square [—6,6] x [—6,6] initially filled 
with constant density and pressure values p — 10^^ and 
p = 5 ■ 10~^ and threaded by a constant horizontal mag- 
netic field with strength Bq — 1. A cylinder with radius 
r = 0.8 and centered at the origin delimits a higher den- 
sity and pressure region where p = 10~^, p = 1. The 
ideal EoS with with F = 4/3 is used. Our choice of 
parameters results in a low /3 plasma {2p/B^ — I0~^) 
and corresponds to the str ongly magnetized cylin drical 

We 
le one 



explosion discussed by Beckwith fc Stone] (2011 

point out th at this con figuration ditt'ers irom t 

discus sed in Komissarov^ ( ^1999) and ^Mignone fc Bodo 
(]2006 1 by having the ambient pressure ten times larger 



Jhis allows to run the computation without having to 
introduce any specific change or modification in the al- 
gorithm, such as shock fiattening or redefinition of the 
total energy. We adopt a base grid of 48 zones in each 
direction with outflow boundary conditions holding on all 
sides. Integratio n proceeds until t = 4 u sing the HLLC 
Riemann solver (Mignone fc Bodo 2006) using 5 levels 
of refinement, reaching an effective resolution of 1536^ 
zones. Total energy triggers refinement with a threshold 
of Xrcf = 0.025. 

The over-pressurized region. Fig. |28[ sets an 
anisotropic blast wave delimited by a weak fast for- 
ward shock propagating almost radially. The explosion 
is strongly confined along the x direction by the mag- 
netic field where plasma is accelerated to 7niax ~ 1.76. 
The inner structure is delimited by an elongated struc- 
ture enclosed by a slow shock adjacent to a contact dis- 
continuity. The two fronts blend together as the prop- 
agation becomes perpendicular to the field line. Since 
the problem is purely two dimensional, rotational dis- 
continuties are absent. The block distribution, shown in 
the lower panel of Fig. [28] shows that refinement clusters 
around the outer discontinuous wave as well as the mul- 
tiple fronts delimiting the extended inner region. The 
total number of blocks is 2258 with the finest level giv- 
ing a relative contribution of '--^ 0.55 and a correspond- 
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Figure 26. The relativistic rotor problem at t = 0.4 using 6 levels of refinement on a base grid with 64^ zones. From left to right, 
top to bottom: density, gas pressure, magnetic field, density, Lorentz factor (magnetic field lines are overplotted) , block distribution and 
domain-average of | V ■ B | . 
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Figure 27. Parallel scalings for the 2D relativistic rotor from 
8 to 512 CPUs. The red and green lines (squares and crosses, re- 
spectively) give the speedup factors corresponding to computations 
with maximum block sizes of 16 and 32 zones, respectively. The 
number of blocks on the finest level {I = 6) at the end of inte- 
gration is reported above and below the corresponding curve. The 
dotted black line gives the perfect scaling. 



ing volume filhng factor of ^ 0.12. The numerical so- 
lution retains the expected degree of symmetry and the 
agreement with earlier results demonstrate the code can 
robustly handle strongly magnetized relativistic shocks 
within the GLM-MHD formalism. Results obtained with 
4 processors show that calculation done with AMR is 
faster than the static grid runs by a factor of ~ 4.5. 



6.5. Spherical Blast Wave 

In the next example we consider a three-dimensional 
ex tension of the cyl i ndrica l blast wave problem presented 



Del Zanna et al. ( 2007 ) . The initial condition consists 
of a sphere with radius r = 0.8 centered at the origin 
filled with hot gas having p = 10~^, p = 1 (the ideal EoS 
with r = 4/3 is used). Outside this region, density and 
pressure decrease linearly reaching the constant ambient 
values p = 10"*^, p — 5 ■ 10^^ for r > 1. The plasma 
is at rest everywhere and a constant magnetic field is 
set oblique to the mesh, i.e., B = Bo{ex + e.y)l\f2 with 
i?o ~ 0.1. The computational domain is the Cartesian 
box [—6, 6]^ covered by a base grid of 40 zones in each di- 
rection with outflow boundary conditions holding on all 
sides. Using the HLLD Riemann solver, we follow inte- 
gration until t — \ using 4 levels of refinement (effective 
resolution 640'^) and a CFL number Ca = 0.4. Zones are 
tagged for refinement upon density fluctuations with a 
threshold value of Xrof = 0-15. 

Results shown in Figure [29] reveal an oval-shaped ex- 
plosion enclosed by an outer fast shock propagating al- 
most radially. A strong reverse shock confines a prolate 
spheroidal region where magnetic field has been drained 
and inside which expansion takes place radially. The 
largest Lorentz factor '-^ 6.3 is attained in the direction 
of magnetic field lines close to the ellipsoid poles. 

In order to measure parallel performance on a fixed 
number of blocks, we have integrated the solution ar- 
ray starting from i = 3.5 for a fixed number of steps 
by enforcing zero interface flux. Results are shown in 
Fig. [30] where we compare two sets of "frozen-flux" com- 
putations with maximum block sizes of 20 and 40 mesh 
points corresponding, respectively, to 27661 (22598 on 
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Figure 28. Cylindrical relativistic blast wave at t = 4. The base level grid has 48^ zones and 5 levels of refinement are used. The difli'erent 
maps show proper density (top left, log scale), thermal pressure (top right) magnetic field strength (top right, log scale) density (bottom 
left, log scale) and Lorentz factor (bottom right). Refinement levels 3 and 4 are overplotted on the first plot while magnetic field lines 
overlap in the latter. 



the finest level) and 11229 (9462 on the finest level) total 
number of blocks. In the former case, the larger number 
of patches allows to reach a better efficiency (more than 
0.75) whereas the latter case performs worse for increas- 
ing number of processors. Even if the finest level has a 
number of blocks still larger than the maximum number 
of CPUs employed, many boxes are much smaller than 
the maximum block size allowed. The ideal workload per 
processor (i.e. the number of cells of the finest level di- 
vided by the number of processors) is comparable to the 
maximum block size in the 1024 CPUs run and becomes 
smaller in the 2048 CPUs run. In this latter case, the 
workload of the processors integrating the larger blocks is 
therefore larger than the ideal one and most of the CPUs 
must wait for these processors to complete the integra- 
tion. This could explain the poor parallel performance 
of this case. 

6.6. Kelvin- Helmholtz flow 
In the next e xample fBucciantini & Del Zanna 2006 



Mignone et al.| 2009) we consider a two-dimensional pla- 
nar domain witiT a; € [0, 1], ?/ S [—1,1] initially filled with 
a (relativistically) hot gas with constant density p = 1 
and pressure p = 20. An initially perturbed shear veloc- 
ity profile of the form 



y 

Vq tanli — , e sin(27ra;) exp 



,0 



r 

is set across the domain, where Vq = 1/4 is the nominal 
flow velocity, e — Uq/IOO is the amplitude of the pertur- 
bation while a ^ 0.01, /3 = 0.1. The TM equation of 
state, Eq. (Ill, is used during the evolution to recover 



the gas enthalpy from density and pressure. The mag- 
netic field is initially uniform and prescribed in terms of 
the components parallel and perpendicular to the plane 
of integration; 



{B.„By,B,) = (^2CTpoip,0,y2^ 



rP 



(74) 



where Cpoi = 0.01, txtor = 1- We perform the integration 
on a base grid of 32 x 64 cells with 6 additional levels of 
refinement, reaching an effective resolution of 2048 x 4096 
zones. The Courant number is set to C„ — 0.8 and re- 
finement is triggered whenever Eq. 
total energy density exceeds 



( 50 1 computed with 
Open bound- 



ary conditions hold at the lower and upper y boundaries 
while peri odicity is imposed in the x direction. 

Fig. [3l] shows density and magnetic field distributions 
at t = 5 which approximately marks the transition from 
the linear phase to the nonline ar evolution (see the dis- 
cussion in Mignone et al.|[2"009[ ). Both density and mag- 
netic field are organized into filaments that wrap around 
forming a number of vortices arranged symmetrically 
with respect to the central point. Refined levels con- 
centrate mainly around the location of the fluid interface 
and accurately follow the formation of the central vortex 
and the more elongated adjacent ones. The computation 
preserves the initial point-symmetry even if the grid gen- 
eration algorithm does not necessarily do so. The lower 
resolution outside this region contributes to damp acous- 
tic waves as they travel towards the outer boundaries and 
eventually reduces the amount of spurious reflections. 

In Fig. [32] we compare the parallel speedup between a 
number of adaptive grid calculations using either 6 lev- 
els (grid jump of 2) or 3 levels (grid jump of 4) and the 
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Figure 29. Density slice-cut for the relativistic magnetized blast wave in three dimensions at t = 4. Four levels of refinement are used to 
achieve an effective resolution of 640'^. The box in the lower half-semispace emphasize jump ratios across levels. Oblique magnetic ficldlines 
are overplotted. 



equivalent uniform grid run with 2048 x 4096 zones. Gen- 
erally, the AMR approach is « 40 faster than the fixed 
grid run. In addition, the 3-level computation seems to 
perform somewhat better than the 6-level run although 
the number of blocks on the finest levels is essentially the 
same at the end of computation (1085 and 1021, respec- 
tively). Parallel scaling sensibly deteriorates when the 
number of blocks per processor becomes less than 3 or 4, 
in accordance with the results of previous tests. 

6.7. Three dimensional Shock-Cloud Interaction 

In the next example we consider a three-dimensional 
extension of the planar relativistic shock-cloud in- 



teraction originally presented in Mignone fc Bodo] 
(20061. The initial condition consists of a high den- 
suy \p = 10) spherical clump (radius 0.15) cen- 
tered at (0.8,0,0) adjacent to a shock wave located at 
X — with downstream and upstream values given by 
{p,v,^,B^,p) = (42.5942, 0, -2.12971, 127.9483) for x < 
and {p,v^,B^,p) = (1, -^/a99, 0.5, lO'^) for x > 0. 
The remaining quantities are set to zero. The computa- 
tional domain is the box a; S [0, 2], y, z S [— i, i]. At 
the base level we fix the resolution to 64 x 32 x 32 zones 
and employ 4 levels of refinement equivalent to an effec- 
tive resolution of 1024 x 512 x 512. The refinement cri- 
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Figure 30. Parallel scalings for the three-dimensional relativistic 
blast wave problem from 32 to 2048 processors. The red and green 
lines (squares and crosses, respectively) give the speedup factors 
corresponding to " frozen- flux" computations with maximum block 
sizes of 20 and 40 zones, respectively. The number of blocks on the 
finest level, reported above and below each line, stays constant in 
time. Ideal scaling is given by the dotted black line. 



terion (501 uses the conserved density for zone tagging 
with a tEreshold Xr = 0-2. The equatio ns of relativistic 
MH D are solved using the TM EoS (Mig none fc McKin- 
ney 2007), the HLLD Riemann solver and linear recon- 
struction with the harmonic limiter (41 ). Shock adaptive 
hybrid integration provides additional numerical dissipa- 
tion in proximity of strong shock waves by switching the 
integration scheme to the HLL Riemann solver and the 
MinMod limiter (Eg . [35]) , according to the strategy de- 
scribed in Appendix |b[ Eor efficiency purposes, we take 
advantage of the symmetry across the xy and xz planes 
to red uce the computation to one quadrant only. 

Fig. 33 shows a three-dimensional rendering of density 
and magnetic pressure through a cross-sectional view of 
the xy and xz planes at i = 1. The collision generates a 
fast, forward bow shock propagating ahead and a back- 
ward reverse shock transmitted back into the cloud. The 
cloud becomes gradually wrapped by the incident shock 
into a mushroom-shaped shell reaching a large compres- 
sion (pmax ~ 121, B^/2 fa 157). Grid refinement con- 
centrates mainly around the incident and forward shocks 
and the tangential discontinuities bordering the edge of 
the cloud. 

6.8. Propagation of Three-dimensional Relativistic 
Magnetized Jet 

As a final application, we discuss the three-dimensional 
propagation of a relativistic magnetized jet carrying an 
initially purely azimuthal magnetic field. Indeed, the 
presence of a substantial toroidal component of the field 
is commonly invoked and held responsible for the accel- 
eration and collimation of jets from active galactic nu- 
clei. If on one hand cylindrical MHD configuration are 
expected to be unstable to reflection, Kelvin-Helmholtz 
and current driven modes, on the other hand astrophysi- 
cal jets a ppear to be quite stabl e thus posing an unsolved 
issue, see Mignone et al. (20101 and reference therein. 



Following Mignone et al. (20091 we set the box x,y ^ 
[—12.5,12.5], z e [0, 50J as our computational domain, 
initially filled with constant uniform density and pressure 
Pa and Pa- The beam is injected at the lower boundary 
for z < 0, r = y/ x^ + y"^ < 1 with a density pj, a longi- 
tudinal velocity corresponding to a Lorentz factor 7j and 
an azimuthal magnetic field given by 



Ijbmf/o, for r<a, 
Brh ~ i Jjbma/r for a < r < 1 
otherwise , 



(75) 



where is the magnetic field strength in the fluid's 
frame and a = 0.5 is the magnetization radius. The 
pressure profile is found by imposing radial momentum 
balance across the beam giving 



1 - min 1^,1 



(76) 



where pa is the ambient pressure deflned in terms of the 
sonic Mach number M: 



Pa 



r(r - i)Af 2 _ Yv^: 



(77) 



with r = 5/3 although the TM equation of state ( 11 ) is 
used during the evolution. The actual value of bm may 
be found by prescribing the magnetization parameter cr^ 
defined as the ratio between beam-averaged magnetic en- 
ergy density and therm al pressure. The final result yields 
( [Mignone et~al]|2009[ ) : 



^ — 



a?{2a^ - 1 + 41oga) 



(78) 



For the present simulation we adopt 7j = 7, M = 3, 
(70 = 1.3, pj = 1 and pa — 10^ thus corresponding to a 
magnetically dominated, underdense relativistic jet. 

We follow the simulation up io t = 130 (in units of 
the speed of light crossing of the jet radius), starting 
with a base grid of 32 x 32 x 64 with 4 levels of refine- 
ment. The refinement ratio between consecutive levels is 
2 so that, at the finest resolution, we have « 20 points 
per beam radius. Zones belonging to levels 0,1 and 2 
are tagged for r efin ement using the normalized second- 
derivative Eq. (50) of the total energy density with a 
threshold Xr = 0^75. Zones at the finest level grid are 
instead created by using 6^/(7^) in place of the total 
energy density. The rationale for choosing this selective 
rule is to provide the largest resolution on the jet mate- 
rial only, still being able to track the sideway expansion 
of the cocoon at lower resolution. We use the HLLD 



Riemann solver with the harmonic limiter ( 41 ) except 
at strong shocks where we employ the multidimensional 
shock dissipation switch outlined in Appendix [B] The 
Courant number is Ca = 0.3. 
The jet structure. Fig. (34), shows that the flow main- 



tains a highly relativistic central spine running through 
multiple recollimation shocks where jet pinching occurs. 
At the jet head the beam decelerates to sub-relativistic 
velocities favoring the formation of a strongly magne- 
tized termination shock where w 267'^p. Here, the 
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Figure 31. Relativistic Kelvin-Helmholtz instability at t = 5 using 6 additional levels of refinement starting from a base grid of 32 X 64. 
The panel on the left shows density (upper half) and the quantity (B^ + By)^ /B^ (lower half) on the whole computational domain. 

7. SUMMARY 

In this paper, we have presented a cell-centered im- 
plementation of the PLUTO code for multi-dimensional 
adaptive mesh refinement (AMR) computations target- 
ing Newtonian and relativistic magnetized flows. A block 
structured approach has been pursued by taking full ad- 
vantage of the high-level parallel distributed infrastruc- 
ture available in the CHOMBO library. This choice 
provides the necessary level of abstraction in delivering 
AMR functionality to the code allowing, at the same 
time, full compatibility with most of the modular imple- 
mentations already available with the static grid version. 
This eases up the process of adding or replacing physics 
modules as long as they comply with the interface re- 
quirements. 

A novel extension to incorporate diffusion terms such 
as viscosity, resistivity and heat conduction without the 
need for operator splitting, has been illustrated in the 
context of the dimensionally unsplit corner transport up- 
wind (CTU) time stepping scheme. In addition, the 
scheme has also been extended to the realm of relativis- 
tic MHD (RMHD) using a characteristic projection-free, 
MUSCL-Hancock normal predictor step. The integration 



Selected regions are enlarged in the two panels on the right. 

presence of a predominant toroidal magnetic field in- 
duces current driven (CD) kink instability which are re- 
sponsible for the observed jet wi ggling and symmetry 



breaking (seelMignone et al. This test shows that 



PLUTO-CHOMBO can be ettectively used in simulations 
involving relativistic velocities, highly supersonic flows 
and strongly magnetized environments. 

At t = 130 the finest level grid has a volume filling fac- 
tor of « 11 per cent with 10761 blocks while level 3 occu- 
pies « 41 per cent of the total volume with 5311 blocks. 
For this particular application, the benefits offered by 
grid adaptivity are more evident at the beginning of the 
computation when most of the computational zones in 
the domain are unrefined. In this respect, we have found 
that the CPU time increases (for t < 130) with the num- 
ber of steps n approximately as a + bn + cr? -t- dr? where 
a K. 1.46, h K. 0.079, c w 5.03 • lO""* and d w 2.38 • 10"^. 
This suggests that, as long as the filling factor remains 
well below 1, the AMR calculation with the proposed re- 
finement criterion should be at least ^ 2.5 times faster 
than an equivalent computation using static mesh refine- 
ment and considerably larger if compared to a uniform 
grid run at the effective resolution. 
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Figure 32. Parallel speedup up to 512 processors for the 
two-dimensional relativistic Kelvin-Helmholtz instability problem. 
Speedup is computed as the ratio /T^^^^ between the ex- 
ecution time of the static uniform grid run at the same effective 
resolution (2048 X 4096) on 8 processors and the running time mea- 
sured with A^CPU processors. Red lines corresponds to the 3-level 
run using max grid sizes of 16 (squares) and 32 (crosses). Green 
lines corresponds the the 6-level compuations. 



scheme retains second-order spatial and temporal accu- 
racy requiring 6 Riemann solvers per cell per step. Inter- 
face states are calculated using the piecewise parabolic 
method (PPM) integration although alternatives based 
on weighted essentially non-oscillatory (WENO) or lin- 
ear slope-limited schemes are also available. The pro- 
posed cell-centered version of PLUTO-CHOMBO en- 
forces the divergence-free condition by augmenting the 
system of equati ons with a generalized Lagrange multi- 
plier (GLM, see |Ded ner et al. 2002 ), in the imp lemen- 
tation outlined by Mignone & Tzeferacos (2010). The 
choice of the GLM-MHD formalism has proven to be a 
convenient starting point in porting a significant frac- 
tion of the static grid code implementations to the AMR 
framework. Although future extensions will also con- 
sider other strategies to enforce the V • B = condition, 



the proposed formulation offers substantial ease of im 
plementation and a viable robust alternative to a con 
strained transport approach which typically requires ad 
ditional care for prolongation and restric tion operations 
at fine-coarse interfaces (see, for instance, Balsara [2001 
Balsara|[2004{ |Cunningham et al._ 20091 



A suite of several test problems including standard nu- 
merical benchmarks and astrophysical applications for 
MHD and RMHD has been selected to assess the effi- 
ciency of PLUTO-CHOMBO in resolving complex flow 
patterns viable through an adaptive grid approach in 1, 
2 and 3 dimensions. Examples include simple shock-tube 
problems, resistive sheets, radiative and thermally con- 
ducting flows, fluid instabilities and blast wave explosions 
in highly magnetized environments. The computational 
saving offered by an adaptive grid computation has been 
shown, for many of the selected test, to be several times 
or even order of magnitudes faster than a static grid in- 
tegration. Parallel performance, evaluated for some of 
the proposed tests, suggests good scalability properties 
for two- and three-dimensional problems as long as the 
number of grids per processor is larger than a factor be- 
tween 2 and 3. 

Both the static and AMR version of PLUTO are 
distributed as a single software package publicly avail- 
able at http://plutocode.ph.unito.it/ while the 
CHOMBO library can be freely downloaded from 
https://seesar.lbl.gov/anag/chombo/. The AMR 
version of PLUTO has been designed to retain salient 



features ch aracterizing the static grid version (Mignone 
et al. 2007) such as modularity - the possibility to eas- 
ily combine different numerical schemes to treat different 
physics-, portability and user- friendliness. While stable 
code releases along with extensive documentation and 
benchmarks are made available on the WEB, the code is 
being actively developed and future development will ad- 
dress additional new physical aspects as well as improved 
numerical algorithms. 

This work has been supported by the PRIN-INAF 
2009 grant. We acknowledge the CINECA Awards N. 
HP10CJ1J54 and N. HPIOBHHHEJ, 2010 under ISCRA 
initiative for the availability of high performance comput- 
ing resources and support. Intensive parallel computa- 
tions were performed on the 4.7 GHz IBM Power 6 p575 
cluster running AIX 6 or the 0.85 GHz IBM BlueGene/P 
cluster running CNL. 



APPENDIX 

DISCRETIZATION OF THE THERMAL CONDUCTION FLUX 



The discretization of the heat conduction flux. 



Eq. (§, 



reflects t he mix ed parabolic/hyperbolic mathematical nature 
of the underlying differential operators, as anticipateH' in Section 2.1.1 From the diffusion limit, |Fciass|/9 — >■ 0, we 



compute Fciass 
x-faces, 



at cell interfaces using second-order accurate central difference expressions for VT, e.g., at constant 



dT 
dx 



Ax 



dT 
dy 



T, 



4Ay 



(Al) 



and similarly at constant y- and z- faces. In the purely hyperbolic limit, |Fciass|/9 — ^ oo, we see that Fc — ^ qt 
where t = Fdass /| Fciass | is a unit vector in the direction of the heat flux. To gain more insights, we re-write the 
one-dimensional energy equation in this limit by keeping only pressure-related terms, that is, 

^ ^ P ^ ^ (50pcLt.)=O, (A2) 



dt 



1- 



1 



dx 
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Figure 33. Vertical cut in the xz plane showing the interaction of a magnetized relativistic shock with a density cloud at t = 1. Density 
and magnetic pressure B^/2 are shown in the left and right panels, respectively. For each quantity, a volumetric rendering is shown in 
the upper half whereas an intensity color map in the middle vertical and horizontal planes displays below. The base grid corresponds to 
64 X 32^ and 4 levels of refinement are employed (effective resolution 1024 X 512^). 



where Cjso — \fvTp the isothermal speed of sound and = • t. Equation (A2) is a nonhnear advection equation 
of the form dtu + dx] = where u = p/{'^ — 1) and / = —qt^ is the hype rbolic saturated flux. A stable discretization 
is therefore provided by adopting an upwind scheme ( Balsara et al.|[2008 ) such as 



= \ 



50 3/2 
Pj 

50 3/2 

Pr 



if < , 

otherwise . 



(A3) 



In the previous equation Pi+i = {pL + Pb)/'2, Pl and pn are computed from the left/right input states for the 

Riemann problem. Finally, we define the thermal conduction flux at constant x-faces (for instance) as the harmonic 
mean between the two regimes, 

iFciassL+i + q., 



(ea; 



(A4) 



where g^+i is computed from Eq. (A3). 



MULTIDIMENSIONAL SHOCK DETECTION AND ADAPTIVE HYBRID-INTEGRATION 

In regions of strong shocks or gradients spurious numerical oscillations may arise and quickly lead to the occurrence 
of unphysical states characterized, for example, by negative pressures, energies, densities or, in the case of relativistic 
flows, superluminal speeds. Such episodes are usually limited to very few grid zones and are originated by either 
insufficient numerical diffusion or lack of a physical solution of the Riemann problem. 

To circumvent this potential hitch, PLUTO provides a safeguard built-in mechanism that i) flag zones that may 
potentially reside inside a strong shock wave and ii) introduces additional numerical dissipation by locally replacing the 
integration scheme with a more diffusive one. For these reasons, most of the interpolation schemes and Riemann solvers 
available with PLUTO embody a hybrid selective mechanism that allows to switch, if required, to a more dissipative 
choice. This process is controlled by a 3D array of integers where each element represents a set of flags that can be 
individually turned on or off by simple bitwise operations. Zones are flagged according to a flexible shock-detection 
criterion, 

|A.g| , \Ayq\ \A,q\ 



min(qi+i,qi_i) mm{qj+i, q^^i) min (q^+i, g^-i) 



V-v<0, 



(Bl) 



where is a standard central difference operator, q is thermal (default) or magnetic pressure, v is the velocity and eg 
is a free adjustable parameter (default is 5). The first condition detects zones within a strong gradient while the second 
one is switched on where compressive motion takes place. When both conditions are met, we flag the zone {i, j, k} to 
be updated with the HLL Riemann solver. At the same time we also flag all neighboring zones whose interpolation 
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Figure 34. Three dimensional visualization of the relativistic magnetically dominated jet at t = 130 using 4 levels of refinement. In the 
top panels we show the distribution of the Lorentz factor (left) and magnetic to kinetic energy ratio |Bp/(7^p) (right). In the bottom 
panels we show the thermal pressure distribution (left) and a slice cut of density together with the grid structure (right). 
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stencil includes the shocked zone to be interpolated using slope-limited reconstruction with the MinMod limiter. We 
note that both flags may be selected and combined independently and, in any case, preserve the second-order accuracy 
of the scheme. 
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